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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/sedini.F90 @ 11504

Last change on this file since 11504 was 11504, checked in by davestorkey, 4 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Strip out all references to nn_dttrc
and the trcsub.F90 module. Notes:

  1. This version of the code currently breaks the GYRE_PISCES test in SETTE.
  2. With the removal of this option, TOP should use the OCE time index variables, eg. Nbb_trc -> Nbb, nittrc000 -> nit0000 etc.
  • Property svn:keywords set to Id
File size: 30.0 KB
Line 
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
12   USE sed_oce
13   USE sedarr
14   USE sedadv
15   USE trcdmp_sed
16   USE trcdta
17   USE iom
18   USE lib_mpp         ! distribued memory computing library
19
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! Module variables
25   REAL(wp)    ::  &
26      sedzmin = 0.3    ,  &  !: Minimum vertical spacing
27      sedhmax = 10.0   ,  &  !: Maximum depth of the sediment
28      sedkth  = 5.0    ,  &  !: Default parameters
29      sedacr  = 3.0          !: Default parameters
30     
31   REAL(wp)    ::  &
32      porsurf =  0.95  ,  &  !: Porosity at the surface
33      porinf  =  0.75  ,  &  !: Porosity at infinite depth
34      rhox    =  2.0         !: Vertical length scale of porosity variation
35
36   REAL(wp)    ::  &
37      rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1]
38      dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s]
39
40   REAL(wp), PUBLIC    ::  &
41      redO2    =  172.  ,  &  !: Redfield coef for Oxygen
42      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate
43      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate
44      redC     =  122.  ,  &  !: Redfield coef for Carbon
45      redfep   =  0.175 ,  &  !: Ratio for iron bound phosphorus
46      rcorgl   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]   
47      rcorgs   =   0.5  ,  &  !: reactivity of the semi-labile component
48      rcorgr   =   1E-4 ,  &  !: reactivity of the refractory component
49      rcnh4    =   10E6 ,  &  !: reactivity for O2/NH4 [l.mol-1.an-1] 
50      rch2s    =   1.E5 ,  &  !: reactivity for O2/ODU [l.mol-1.an-1]
51      rcfe2    =   5.E8 ,  &  !: reactivity for O2/Fe2+ [l.mol-1.an-1]
52      rcfeh2s  =   1.E4 ,  &  !: Reactivity for FEOH/H2S [l.mol-1.an-1]
53      rcfes    =   1.E5 ,  &  !: Reactivity for FE2+/H2S [l.mol-1.an-1]
54      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1]
55      xksedo2  =   5E-6 ,  &  !: half-sturation constant for oxic remin.
56      xksedno3 =   5E-6 ,  &  !: half-saturation constant for denitrification
57      xksedfeo =   0.6  ,  & !: half-saturation constant for iron remin
58      xksedso4 =   2E-3       !: half-saturation constant for SO4 remin
59
60   REAL(wp)    ::  &
61      rccal   = 1000.,      & !: reactivity for calcite         [l.mol-1.an-1]
62      rcligc  = 1.E-4         !: L/C ratio in POC
63
64   REAL(wp), PUBLIC    ::  dbiot   = 15. , &  !: coefficient for bioturbation    [cm**2.(n-1)]
65      dbtbzsc =  10.0  ,    &  !: Vertical scale of variation. If no variation, mixed layer in the sed [cm]
66      xirrzsc = 2.0            !: Vertical scale of irrigation variation.
67   REAL(wp)    ::  &
68      ryear = 365. * 24. * 3600. !:  1 year converted in second
69
70   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff1
71   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 /
72
73   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2
74   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 /
75
76
77
78   !! *  Routine accessibility
79   PUBLIC sed_init          ! routine called by opa.F90
80
81   !! $Id$
82CONTAINS
83
84
85   SUBROUTINE sed_init
86      !!----------------------------------------------------------------------
87      !!                   ***  ROUTINE sed_init  ***
88      !!
89      !! ** Purpose :  Initialization of sediment module
90      !!               - Reading namelist
91      !!               - Read the deepest water layer thickness
92      !!                 ( using as mask ) in Netcdf file
93      !!               - Convert unity if necessary
94      !!               - sets initial sediment composition
95      !!                 ( only clay or reading restart file )
96      !!               - sets sediment grid, porosity and others constants
97      !!
98      !!   History :
99      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
100      !!        !  06-07  (C. Ethe)  Re-organization
101      !!----------------------------------------------------------------------
102      INTEGER :: ji, jj, ikt, ierr
103      !!----------------------------------------------------------------------
104
105
106      ! Reading namelist.sed variables
107      !---------------------------------------
108
109      CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
110
111      IF (lwp) THEN
112         WRITE(numsed,*)
113         WRITE(numsed,*) '                 PISCES framework'
114         WRITE(numsed,*) '                 SEDIMENT model'
115         WRITE(numsed,*) '                version 3.0  (2018) '
116         WRITE(numsed,*)
117         WRITE(numsed,*)
118      ENDIF
119
120      IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
121      IF(lwp) WRITE(numsed,*) ' '
122
123      ! Read sediment Namelist
124      !-------------------------
125      CALL sed_init_nam
126
127      ! Allocate SEDIMENT arrays
128      ierr =        sed_alloc()
129      ierr = ierr + sed_oce_alloc()
130      ierr = ierr + sed_adv_alloc() 
131      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' )
132
133      ! Determination of sediments number of points and allocate global variables
134      epkbot(:,:) = 0.
135      DO jj = 1, jpj
136         DO ji = 1, jpi
137            ikt = mbkt(ji,jj) 
138            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt)
139            gdepbot(ji,jj) = gdepw_0(ji,jj,ikt)
140         ENDDO
141      ENDDO
142
143      ! computation of total number of ocean points
144      !--------------------------------------------
145      sedmask = 0.
146      IF ( COUNT( epkbot(:,:) > 0. ) == 0 ) THEN
147          sedmask = 0.
148      ELSE
149          sedmask = 1.
150      ENDIF
151      jpoce  = MAX( COUNT( epkbot(:,:) > 0. ) , 1 )
152
153      ! Allocate memory size of global variables
154      ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) )
155      ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ;  ALLOCATE( rainrm_dta(jpoce,jpsol) )
156      ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) ) 
157      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) )
158      ALLOCATE( zkbot(jpoce) )               ;  ALLOCATE( db(jpoce,jpksed) )
159      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) ) 
160      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) )
161      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) )
162      ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) ) 
163      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) )
164      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) )
165      ALLOCATE( sedligand(jpoce, jpksed) )
166
167      ! Initialization of global variables
168      pwcp  (:,:,:) = 0.   ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0. 
169      solcp (:,:,:) = 0.   ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0.
170      rainrm(:,:  ) = 0.   ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0. 
171      dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0.
172      temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0.
173      press (:    ) = 0.   ;  densSW (:   ) = 0.  ; db        (:,:) = 0. 
174      hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0. 
175      dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0. 
176      fecratio(:)   = 1E-5 
177      sedligand(:,:) = 0.6E-9
178
179      ! Chemical variables     
180      ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) )     
181      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )   
182      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )   
183      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) ) 
184      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) )
185      ALLOCATE( fluorids(jpoce) )
186
187      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0.
188      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0.
189      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0.
190      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0.
191      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0.
192
193      ! Mass balance calculation 
194      ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ;  ALLOCATE( rloss(jpoce, jpsol) )
195      ALLOCATE( tokbot (jpoce, jpwat) ) 
196
197      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0. 
198
199      ! Initialization of sediment geometry
200      !------------------------------------
201      CALL sed_init_geom
202
203      ! Offline specific mode
204      ! ---------------------
205      ln_sediment_offline = .FALSE.
206
207#if defined key_sed_off
208      ln_sediment_offline = .TRUE.
209      IF (lwp) write(numsed,*) 'Sediment module is run in offline mode'
210      IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time'
211      IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false'
212      IF (lwp) write(numsed,*) '--------------------------------------------'
213      ln_sed_2way = .FALSE.
214#endif
215      ! Initialisation of tracer damping
216      ! --------------------------------
217      IF (ln_sediment_offline) THEN
218         CALL trc_dta_ini(jptra)
219         CALL trc_dmp_sed_ini
220      ENDIF
221
222   END SUBROUTINE sed_init
223
224   SUBROUTINE sed_init_geom
225      !!----------------------------------------------------------------------
226      !!                   ***  ROUTINE sed_init_geom  ***
227      !!
228      !! ** Purpose :  Initialization of sediment geometry
229      !!               - Read the deepest water layer thickness
230      !!                 ( using as mask ) in Netcdf file
231      !!               - sets sediment grid, porosity and molecular weight
232      !!                 and others constants
233      !!
234      !!   History :
235      !!        !  06-07  (C. Ethe)  Original
236      !!----------------------------------------------------------------------
237      !! * Modules used
238      !! * local declarations
239      INTEGER  :: ji, jj, jk, jn
240      REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw
241      REAL(wp) :: ztmp1, ztmp2
242      !----------------------------------------------------------
243
244      IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
245      IF(lwp) WRITE(numsed,*) ' '
246
247      ! Computation of 1D array of sediments points
248      indoce = 0
249      DO jj = 1, jpj
250         DO ji = 1, jpi
251            IF (  epkbot(ji,jj) > 0. ) THEN
252               indoce          = indoce + 1
253               iarroce(indoce) = (jj - 1) * jpi + ji
254            ENDIF
255         END DO
256      END DO
257
258      IF ( indoce .EQ. 0 ) THEN
259         indoce = 1
260         iarroce(indoce) = 1
261      ENDIF
262
263      IF( indoce .NE. jpoce ) THEN
264         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' )
265      ELSE
266         IF (lwp) WRITE(numsed,*) ' '
267         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
268         IF (lwp) WRITE(numsed,*) ' '
269      ENDIF
270
271      ! initialization of dzkbot in [cm]
272      !------------------------------------------------   
273      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
274      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
275      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) )
276
277      ! Geometry and  constants
278      ! sediment layer thickness [cm]
279      ! (1st layer= diffusive layer = pur water)
280      !------------------------------------------
281      za1  = (  sedzmin - sedhmax / FLOAT(jpksed-1)  )                                                      &
282         & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * (  LOG( COSH( (jpksed - sedkth) / sedacr) )      &
283         &                                                   - LOG( COSH( ( 1  - sedkth) / sedacr) )  )  )
284      za0  = sedzmin - za1 * TANH( (1-sedkth) / sedacr )
285      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  )
286
287      profsedw(1) = 0.0
288      profsed(1) = -dz(1) / 2.
289      DO jk = 2, jpksed
290         zw = REAL( jk , wp )
291         zt = REAL( jk , wp ) - 0.5_wp
292         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  ) 
293         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  )
294      END DO
295
296      dz(1) = 0.1
297      DO jk = 2, jpksed
298         dz(jk) = profsedw(jk) - profsedw(jk-1)
299      END DO
300
301      DO jk = 1, jpksed
302         DO ji = 1, jpoce
303            dz3d(ji,jk) = dz(jk)
304         END DO
305      ENDDO
306
307      !  Porosity profile [0]
308      !---------------------
309      por(1) = 1.0
310      DO jk = 2, jpksed
311         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) )
312      END DO
313 
314      ! inverse of  Porosity profile
315      !-----------------------------
316      por1(:) = 1. - por(:)
317
318      ! Volumes of pore water and solid fractions (vector and array)
319      !     WARNING : volw(1) and vols(1) are sublayer volums
320      volw(:) = dz(:) * por(:)
321      vols(:) = dz(:) * por1(:)
322
323      ! temporary new value for dz3d(:,1)
324      dz3d(1:jpoce,1) = dzkbot(1:jpoce)
325
326      ! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
327      DO jk = 1, jpksed
328         volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
329         vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
330      ENDDO
331
332      ! Back to the old sublayer vlaue for dz3d(:,1)
333      dz3d(1:jpoce,1) = dz(1)
334
335      !---------------------------------------------
336      ! Molecular weight [g/mol] for solid species
337      !---------------------------------------------
338
339      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
340      !---------------------------------------
341      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 
342
343      !  clay
344      !  some kind of Illit (according to Pape)
345      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]
346      !--------------------------------------------------------------------
347      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
348         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    &
349         &              0.59 * 27. + 10. * 16.
350
351      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0)
352
353      mol_wgt(jsfes)  = 55.0 + 32.0
354
355      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
356      ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
357      ! So C(122)H(355)O(120)N(16)P(1)
358      !------------------------------------------------------------
359      mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
360         &                16. * 14. + 31. ) / 122.
361      mol_wgt(jspos) = mol_wgt(jspoc)
362      mol_wgt(jspor) = mol_wgt(jspoc)
363
364      ! CaCO3
365      !---------
366      mol_wgt(jscal) = 40. + 12. + 3. * 16.
367
368      ! Density of solid material in sediment [g/cm**3]
369      !------------------------------------------------
370      denssol = 2.6
371
372      ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
373      !--------------------------------------------------------------------
374!      DO jn = 1, jpsol
375!         DO jk = 1, jpksed
376!            DO ji = 1, jpoce
377!               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) )
378!            END DO
379!         END DO
380!      END DO
381
382      ! Accumulation rate from Burwicz et al. (2011). This is used to
383      ! compute the flux of clays and minerals
384      ! --------------------------------------------------------------
385      DO ji = 1, jpoce
386          ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 )
387          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 )
388          wacc(ji) = ztmp1 + ztmp2
389      END DO
390
391
392      ! Initialization of time step as function of porosity [cm**2/s]
393      !------------------------------------------------------------------
394   END SUBROUTINE sed_init_geom
395
396   SUBROUTINE sed_init_nam
397      !!----------------------------------------------------------------------
398      !!                   ***  ROUTINE sed_init_nam  ***
399      !!
400      !! ** Purpose :  Initialization of sediment geometry
401      !!               - Reading namelist and defines constants variables
402      !!
403      !!   History :
404      !!        !  06-07  (C. Ethe)  Original
405      !!----------------------------------------------------------------------
406
407      INTEGER ::   numnamsed_ref = -1           !! Logical units for namelist sediment
408      INTEGER ::   numnamsed_cfg = -1           !! Logical units for namelist sediment
409      INTEGER :: ios                 ! Local integer output status for namelist read
410      CHARACTER(LEN=20)   ::   clname
411
412      TYPE PSED
413         CHARACTER(len = 20)  :: snamesed   !: short name
414         CHARACTER(len = 80 ) :: lnamesed   !: long name
415         CHARACTER(len = 20 ) :: unitsed    !: unit
416      END TYPE PSED
417
418      TYPE(PSED) , DIMENSION(jpsol     ) :: sedsol
419      TYPE(PSED) , DIMENSION(jpwat     ) :: sedwat
420      TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
421      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
422
423      NAMELIST/nam_run/nrseddt,ln_sed_2way
424      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox
425      NAMELIST/nam_trased/sedsol, sedwat
426      NAMELIST/nam_diased/seddiag3d, seddiag2d
427      NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc
428      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  &
429         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso,  &
430         &             xksedo2, xksedno3, xksedfeo, xksedso4
431      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc
432      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out
433
434      INTEGER :: ji, jn, jn1
435      !-------------------------------------------------------
436
437      IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists '
438      IF(lwp) WRITE(numsed,*) ' '
439
440      ! ryear = 1 year converted in second
441      !------------------------------------
442      IF (lwp) THEN
443         WRITE(numsed,*) ' '
444         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
445         WRITE(numsed,*) ' '     
446      ENDIF
447
448      ! Reading namelist.sed variables
449      !---------------------------------
450      clname = 'namelist_sediment'
451      IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist'
452      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~'
453      CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
454      CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
455
456      nitsed000 = nittrc000
457      nitsedend = nitend
458      ! Namelist nam_run
459      REWIND( numnamsed_ref )              ! Namelist nam_run in reference namelist : Pisces variables
460      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901)
461901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist', lwp )
462
463      REWIND( numnamsed_cfg )              ! Namelist nam_run in reference namelist : Pisces variables
464      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902)
465902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist', lwp )
466
467      IF (lwp) THEN
468         WRITE(numsed,*) ' namelist nam_run'
469         WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt
470         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way
471      ENDIF
472
473      IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' )
474
475      REWIND( numnamsed_ref )              ! Namelist nam_geom in reference namelist : Pisces variables
476      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
477903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist', lwp )
478
479      REWIND( numnamsed_cfg )              ! Namelist nam_geom in reference namelist : Pisces variables
480      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904)
481904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist', lwp )
482
483      IF (lwp) THEN
484         WRITE(numsed,*) ' namelist nam_geom'
485         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed
486         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin
487         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax
488         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth
489         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr
490         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf
491         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf
492         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox
493      ENDIF
494
495      jpksedm1  = jpksed - 1
496      dtsed = r2dttrc
497
498      REWIND( numnamsed_ref )              ! Namelist nam_trased in reference namelist : Pisces variables
499      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905)
500905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist', lwp )
501
502      REWIND( numnamsed_cfg )              ! Namelist nam_trased in reference namelist : Pisces variables
503      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906)
504906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist', lwp )
505
506      DO jn = 1, jpsol
507         sedtrcd(jn) = sedsol(jn)%snamesed
508         sedtrcl(jn) = sedsol(jn)%lnamesed
509         sedtrcu(jn) = sedsol(jn)%unitsed
510      END DO
511
512      DO jn = 1, jpwat
513         jn1 = jn + jpsol
514         sedtrcd(jn1) = sedwat(jn)%snamesed
515         sedtrcl(jn1) = sedwat(jn)%lnamesed
516         sedtrcu(jn1) = sedwat(jn)%unitsed
517      END DO
518
519      IF (lwp) THEN
520         WRITE(numsed,*) ' namelist nam_trased'
521         WRITE(numsed,*) ' '
522         DO jn = 1, jptrased
523            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
524            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
525            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
526            WRITE(numsed,*) ' '
527         END DO
528         WRITE(numsed,*) ' '
529      ENDIF
530
531      REWIND( numnamsed_ref )              ! Namelist nam_diased in reference namelist : Pisces variables
532      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907)
533907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist', lwp )
534
535      REWIND( numnamsed_cfg )              ! Namelist nam_diased in reference namelist : Pisces variables
536      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908)
537908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist', lwp )
538     
539      DO jn = 1, jpdia3dsed
540         seddia3d(jn) = seddiag3d(jn)%snamesed
541         seddia3l(jn) = seddiag3d(jn)%lnamesed
542         seddia3u(jn) = seddiag3d(jn)%unitsed
543      END DO
544
545      DO jn = 1, jpdia2dsed
546         seddia2d(jn) = seddiag2d(jn)%snamesed
547         seddia2l(jn) = seddiag2d(jn)%lnamesed
548         seddia2u(jn) = seddiag2d(jn)%unitsed
549      END DO
550
551      IF (lwp) THEN
552         WRITE(numsed,*) ' namelist nam_diased'
553         WRITE(numsed,*) ' '
554         DO jn = 1, jpdia3dsed
555            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
556            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
557            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
558            WRITE(numsed,*) ' '
559         END DO
560
561         DO jn = 1, jpdia2dsed
562            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
563            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
564            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
565            WRITE(numsed,*) ' '
566         END DO
567
568         WRITE(numsed,*) ' '
569      ENDIF
570
571      ! Inorganic chemistry parameters
572      !----------------------------------
573      REWIND( numnamsed_ref )              ! Namelist nam_inorg in reference namelist : Pisces variables
574      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909)
575909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist', lwp )
576
577      REWIND( numnamsed_cfg )              ! Namelist nam_inorg in reference namelist : Pisces variables
578      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910)
579910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist', lwp )
580
581      IF (lwp) THEN
582         WRITE(numsed,*) ' namelist nam_inorg'
583         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal
584         WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef
585         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal
586         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc
587         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc
588         WRITE(numsed,*) ' '
589      ENDIF
590
591      ! Unity conversion to get saturation conc. psat in [mol.l-1]
592      ! and reactivity rc in  [l.mol-1.s-1]
593      !----------------------------------------------------------
594      reac_sil   = rcopal / ryear     
595      reac_ligc  = rcligc / ryear
596
597      ! Additional parameter linked to POC/O2/No3/Po4
598      !----------------------------------------------
599      REWIND( numnamsed_ref )              ! Namelist nam_poc in reference namelist : Pisces variables
600      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911)
601911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist', lwp )
602
603      REWIND( numnamsed_cfg )              ! Namelist nam_poc in reference namelist : Pisces variables
604      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912)
605912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist', lwp )
606
607      IF (lwp) THEN
608         WRITE(numsed,*) ' namelist nam_poc'
609         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2
610         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3
611         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4
612         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC
613         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep
614         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl
615         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs
616         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr
617         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4
618         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s
619         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2
620         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s
621         WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes
622         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso
623         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2
624         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3
625         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo
626         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4
627         WRITE(numsed,*) ' '
628      ENDIF
629
630
631      so2ut  = redO2    / redC
632      srno3  = redNo3   / redC
633      spo4r  = redPo4   / redC
634      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC
635      ! reactivity rc in  [l.mol-1.s-1]
636      reac_pocl  = rcorgl / ryear
637      reac_pocs  = rcorgs / ryear
638      reac_pocr  = rcorgr / ryear
639      reac_nh4   = rcnh4  / ryear
640      reac_h2s   = rch2s  / ryear
641      reac_fe2   = rcfe2  / ryear
642      reac_feh2s = rcfeh2s/ ryear
643      reac_fes   = rcfes  / ryear
644      reac_feso  = rcfeso / ryear
645
646      ! reactivity rc in  [l.mol-1.s-1]     
647      reac_cal = rccal / ryear
648
649      ! Bioturbation parameter
650      !------------------------
651      REWIND( numnamsed_ref )              ! Namelist nam_btb in reference namelist : Pisces variables
652      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913)
653913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist', lwp )
654
655      REWIND( numnamsed_cfg )              ! Namelist nam_btb in reference namelist : Pisces variables
656      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914)
657914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist', lwp )
658
659      IF (lwp) THEN
660         WRITE(numsed,*) ' namelist nam_btb ' 
661         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot
662         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz
663         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc
664         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4
665         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig
666         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc
667         WRITE(numsed,*) ' '
668      ENDIF
669
670      ! Initial value (t=0) for sediment pore water and solid components
671      !----------------------------------------------------------------
672      REWIND( numnamsed_ref )              ! Namelist nam_rst in reference namelist : Pisces variables
673      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915)
674915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist', lwp )
675
676      REWIND( numnamsed_cfg )              ! Namelist nam_rst in reference namelist : Pisces variables
677      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916)
678916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist', lwp )
679
680      IF (lwp) THEN
681         WRITE(numsed,*) ' namelist  nam_rst ' 
682         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
683         WRITE(numsed,*) ' '
684      ENDIF
685      nn_dtsed = 1
686
687      CLOSE( numnamsed_cfg )
688      CLOSE( numnamsed_ref )
689
690   END SUBROUTINE sed_init_nam
691
692END MODULE sedini
Note: See TracBrowser for help on using the repository browser.