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/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedini.F90 @ 12460

Last change on this file since 12460 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 28.2 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   !! * Substitutions
25#  include "do_loop_substitute.h90"
26   !! Module variables
27   REAL(wp)    ::  &
28      sedzmin = 0.3    ,  &  !: Minimum vertical spacing
29      sedhmax = 10.0   ,  &  !: Maximum depth of the sediment
30      sedkth  = 5.0    ,  &  !: Default parameters
31      sedacr  = 3.0          !: Default parameters
32     
33   REAL(wp)    ::  &
34      porsurf =  0.95  ,  &  !: Porosity at the surface
35      porinf  =  0.75  ,  &  !: Porosity at infinite depth
36      rhox    =  2.0         !: Vertical length scale of porosity variation
37
38   REAL(wp)    ::  &
39      rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1]
40      dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s]
41
42   REAL(wp), PUBLIC    ::  &
43      redO2    =  172.  ,  &  !: Redfield coef for Oxygen
44      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate
45      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate
46      redC     =  122.  ,  &  !: Redfield coef for Carbon
47      redfep   =  0.175 ,  &  !: Ratio for iron bound phosphorus
48      rcorgl   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]   
49      rcorgs   =   0.5  ,  &  !: reactivity of the semi-labile component
50      rcorgr   =   1E-4 ,  &  !: reactivity of the refractory component
51      rcnh4    =   10E6 ,  &  !: reactivity for O2/NH4 [l.mol-1.an-1] 
52      rch2s    =   1.E5 ,  &  !: reactivity for O2/ODU [l.mol-1.an-1]
53      rcfe2    =   5.E8 ,  &  !: reactivity for O2/Fe2+ [l.mol-1.an-1]
54      rcfeh2s  =   1.E4 ,  &  !: Reactivity for FEOH/H2S [l.mol-1.an-1]
55      rcfes    =   1.E5 ,  &  !: Reactivity for FE2+/H2S [l.mol-1.an-1]
56      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1]
57      xksedo2  =   5E-6 ,  &  !: half-sturation constant for oxic remin.
58      xksedno3 =   5E-6 ,  &  !: half-saturation constant for denitrification
59      xksedfeo =   0.6  ,  & !: half-saturation constant for iron remin
60      xksedso4 =   2E-3       !: half-saturation constant for SO4 remin
61
62   REAL(wp)    ::  &
63      rccal   = 1000.,      & !: reactivity for calcite         [l.mol-1.an-1]
64      rcligc  = 1.E-4         !: L/C ratio in POC
65
66   REAL(wp), PUBLIC    ::  dbiot   = 15. , &  !: coefficient for bioturbation    [cm**2.(n-1)]
67      dbtbzsc =  10.0  ,    &  !: Vertical scale of variation. If no variation, mixed layer in the sed [cm]
68      xirrzsc = 2.0            !: Vertical scale of irrigation variation.
69   REAL(wp)    ::  &
70      ryear = 365. * 24. * 3600. !:  1 year converted in second
71
72   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff1
73   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 /
74
75   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2
76   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 /
77
78
79
80   !! *  Routine accessibility
81   PUBLIC sed_init          ! routine called by opa.F90
82
83   !! $Id$
84CONTAINS
85
86
87   SUBROUTINE sed_init
88      !!----------------------------------------------------------------------
89      !!                   ***  ROUTINE sed_init  ***
90      !!
91      !! ** Purpose :  Initialization of sediment module
92      !!               - Reading namelist
93      !!               - Read the deepest water layer thickness
94      !!                 ( using as mask ) in Netcdf file
95      !!               - Convert unity if necessary
96      !!               - sets initial sediment composition
97      !!                 ( only clay or reading restart file )
98      !!               - sets sediment grid, porosity and others constants
99      !!
100      !!   History :
101      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
102      !!        !  06-07  (C. Ethe)  Re-organization
103      !!----------------------------------------------------------------------
104      INTEGER :: ji, jj, ikt, ierr
105      !!----------------------------------------------------------------------
106
107
108      ! Reading namelist.sed variables
109      !---------------------------------------
110
111      CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
112
113      IF (lwp) THEN
114         WRITE(numsed,*)
115         WRITE(numsed,*) '                 PISCES framework'
116         WRITE(numsed,*) '                 SEDIMENT model'
117         WRITE(numsed,*) '                version 3.0  (2018) '
118         WRITE(numsed,*)
119         WRITE(numsed,*)
120      ENDIF
121
122      IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
123      IF(lwp) WRITE(numsed,*) ' '
124
125      ! Read sediment Namelist
126      !-------------------------
127      CALL sed_init_nam
128
129      ! Allocate SEDIMENT arrays
130      ierr =        sed_alloc()
131      ierr = ierr + sed_oce_alloc()
132      ierr = ierr + sed_adv_alloc() 
133      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' )
134
135      ! Determination of sediments number of points and allocate global variables
136      epkbot(:,:) = 0.
137      DO_2D_11_11
138         ikt = mbkt(ji,jj) 
139         IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt)
140         gdepbot(ji,jj) = gdepw_0(ji,jj,ikt)
141      END_2D
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_2D_11_11
250         IF (  epkbot(ji,jj) > 0. ) THEN
251            indoce          = indoce + 1
252            iarroce(indoce) = (jj - 1) * jpi + ji
253         ENDIF
254      END_2D
255
256      IF ( indoce .EQ. 0 ) THEN
257         indoce = 1
258         iarroce(indoce) = 1
259      ENDIF
260
261      IF( indoce .NE. jpoce ) THEN
262         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' )
263      ELSE
264         IF (lwp) WRITE(numsed,*) ' '
265         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
266         IF (lwp) WRITE(numsed,*) ' '
267      ENDIF
268
269      ! initialization of dzkbot in [cm]
270      !------------------------------------------------   
271      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
272      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
273      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) )
274
275      ! Geometry and  constants
276      ! sediment layer thickness [cm]
277      ! (1st layer= diffusive layer = pur water)
278      !------------------------------------------
279      za1  = (  sedzmin - sedhmax / FLOAT(jpksed-1)  )                                                      &
280         & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * (  LOG( COSH( (jpksed - sedkth) / sedacr) )      &
281         &                                                   - LOG( COSH( ( 1  - sedkth) / sedacr) )  )  )
282      za0  = sedzmin - za1 * TANH( (1-sedkth) / sedacr )
283      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  )
284
285      profsedw(1) = 0.0
286      profsed(1) = -dz(1) / 2.
287      DO jk = 2, jpksed
288         zw = REAL( jk , wp )
289         zt = REAL( jk , wp ) - 0.5_wp
290         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  ) 
291         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  )
292      END DO
293
294      dz(1) = 0.1
295      DO jk = 2, jpksed
296         dz(jk) = profsedw(jk) - profsedw(jk-1)
297      END DO
298
299      DO jk = 1, jpksed
300         DO ji = 1, jpoce
301            dz3d(ji,jk) = dz(jk)
302         END DO
303      ENDDO
304
305      !  Porosity profile [0]
306      !---------------------
307      por(1) = 1.0
308      DO jk = 2, jpksed
309         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) )
310      END DO
311 
312      ! inverse of  Porosity profile
313      !-----------------------------
314      por1(:) = 1. - por(:)
315
316      ! Volumes of pore water and solid fractions (vector and array)
317      !     WARNING : volw(1) and vols(1) are sublayer volums
318      volw(:) = dz(:) * por(:)
319      vols(:) = dz(:) * por1(:)
320
321      ! temporary new value for dz3d(:,1)
322      dz3d(1:jpoce,1) = dzkbot(1:jpoce)
323
324      ! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
325      DO jk = 1, jpksed
326         volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
327         vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
328      ENDDO
329
330      ! Back to the old sublayer vlaue for dz3d(:,1)
331      dz3d(1:jpoce,1) = dz(1)
332
333      !---------------------------------------------
334      ! Molecular weight [g/mol] for solid species
335      !---------------------------------------------
336
337      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
338      !---------------------------------------
339      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 
340
341      !  clay
342      !  some kind of Illit (according to Pape)
343      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]
344      !--------------------------------------------------------------------
345      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
346         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    &
347         &              0.59 * 27. + 10. * 16.
348
349      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0)
350
351      mol_wgt(jsfes)  = 55.0 + 32.0
352
353      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
354      ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
355      ! So C(122)H(355)O(120)N(16)P(1)
356      !------------------------------------------------------------
357      mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
358         &                16. * 14. + 31. ) / 122.
359      mol_wgt(jspos) = mol_wgt(jspoc)
360      mol_wgt(jspor) = mol_wgt(jspoc)
361
362      ! CaCO3
363      !---------
364      mol_wgt(jscal) = 40. + 12. + 3. * 16.
365
366      ! Density of solid material in sediment [g/cm**3]
367      !------------------------------------------------
368      denssol = 2.6
369
370      ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
371      !--------------------------------------------------------------------
372!      DO jn = 1, jpsol
373!         DO jk = 1, jpksed
374!            DO ji = 1, jpoce
375!               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) )
376!            END DO
377!         END DO
378!      END DO
379
380      ! Accumulation rate from Burwicz et al. (2011). This is used to
381      ! compute the flux of clays and minerals
382      ! --------------------------------------------------------------
383      DO ji = 1, jpoce
384          ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 )
385          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 )
386          wacc(ji) = ztmp1 + ztmp2
387      END DO
388
389
390      ! Initialization of time step as function of porosity [cm**2/s]
391      !------------------------------------------------------------------
392   END SUBROUTINE sed_init_geom
393
394   SUBROUTINE sed_init_nam
395      !!----------------------------------------------------------------------
396      !!                   ***  ROUTINE sed_init_nam  ***
397      !!
398      !! ** Purpose :  Initialization of sediment geometry
399      !!               - Reading namelist and defines constants variables
400      !!
401      !!   History :
402      !!        !  06-07  (C. Ethe)  Original
403      !!----------------------------------------------------------------------
404
405      CHARACTER(:), ALLOCATABLE ::   numnamsed_ref           !! Character buffer for reference namelist sediment
406      CHARACTER(:), ALLOCATABLE ::   numnamsed_cfg           !! Character buffer for configuration namelist sediment
407      INTEGER :: ios                 ! Local integer output status for namelist read
408      CHARACTER(LEN=20)   ::   clname
409
410      TYPE PSED
411         CHARACTER(len = 20)  :: snamesed   !: short name
412         CHARACTER(len = 80 ) :: lnamesed   !: long name
413         CHARACTER(len = 20 ) :: unitsed    !: unit
414      END TYPE PSED
415
416      TYPE(PSED) , DIMENSION(jpsol     ) :: sedsol
417      TYPE(PSED) , DIMENSION(jpwat     ) :: sedwat
418      TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
419      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
420
421      NAMELIST/nam_run/nrseddt,ln_sed_2way
422      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox
423      NAMELIST/nam_trased/sedsol, sedwat
424      NAMELIST/nam_diased/seddiag3d, seddiag2d
425      NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc
426      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  &
427         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso,  &
428         &             xksedo2, xksedno3, xksedfeo, xksedso4
429      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc
430      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out
431
432      INTEGER :: ji, jn, jn1
433      !-------------------------------------------------------
434
435      IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists '
436      IF(lwp) WRITE(numsed,*) ' '
437
438      ! ryear = 1 year converted in second
439      !------------------------------------
440      IF (lwp) THEN
441         WRITE(numsed,*) ' '
442         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
443         WRITE(numsed,*) ' '     
444      ENDIF
445
446      ! Reading namelist.sed variables
447      !---------------------------------
448      clname = 'namelist_sediment'
449      IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist'
450      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~'
451      CALL load_nml( numnamsed_ref, TRIM( clname )//'_ref', numout, lwm )
452      CALL load_nml( numnamsed_cfg, TRIM( clname )//'_cfg', numout, lwm )
453
454      nitsed000 = nittrc000
455      nitsedend = nitend
456      ! Namelist nam_run
457      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901)
458901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' )
459
460      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902)
461902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' )
462
463      IF (lwp) THEN
464         WRITE(numsed,*) ' namelist nam_run'
465         WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt
466         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way
467      ENDIF
468
469      IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' )
470
471      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
472903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' )
473
474      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904)
475904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' )
476
477      IF (lwp) THEN
478         WRITE(numsed,*) ' namelist nam_geom'
479         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed
480         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin
481         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax
482         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth
483         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr
484         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf
485         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf
486         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox
487      ENDIF
488
489      jpksedm1  = jpksed - 1
490      dtsed = r2dttrc
491
492      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905)
493905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' )
494
495      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906)
496906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' )
497
498      DO jn = 1, jpsol
499         sedtrcd(jn) = sedsol(jn)%snamesed
500         sedtrcl(jn) = sedsol(jn)%lnamesed
501         sedtrcu(jn) = sedsol(jn)%unitsed
502      END DO
503
504      DO jn = 1, jpwat
505         jn1 = jn + jpsol
506         sedtrcd(jn1) = sedwat(jn)%snamesed
507         sedtrcl(jn1) = sedwat(jn)%lnamesed
508         sedtrcu(jn1) = sedwat(jn)%unitsed
509      END DO
510
511      IF (lwp) THEN
512         WRITE(numsed,*) ' namelist nam_trased'
513         WRITE(numsed,*) ' '
514         DO jn = 1, jptrased
515            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
516            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
517            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
518            WRITE(numsed,*) ' '
519         END DO
520         WRITE(numsed,*) ' '
521      ENDIF
522
523      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907)
524907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' )
525
526      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908)
527908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' )
528     
529      DO jn = 1, jpdia3dsed
530         seddia3d(jn) = seddiag3d(jn)%snamesed
531         seddia3l(jn) = seddiag3d(jn)%lnamesed
532         seddia3u(jn) = seddiag3d(jn)%unitsed
533      END DO
534
535      DO jn = 1, jpdia2dsed
536         seddia2d(jn) = seddiag2d(jn)%snamesed
537         seddia2l(jn) = seddiag2d(jn)%lnamesed
538         seddia2u(jn) = seddiag2d(jn)%unitsed
539      END DO
540
541      IF (lwp) THEN
542         WRITE(numsed,*) ' namelist nam_diased'
543         WRITE(numsed,*) ' '
544         DO jn = 1, jpdia3dsed
545            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
546            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
547            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
548            WRITE(numsed,*) ' '
549         END DO
550
551         DO jn = 1, jpdia2dsed
552            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
553            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
554            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
555            WRITE(numsed,*) ' '
556         END DO
557
558         WRITE(numsed,*) ' '
559      ENDIF
560
561      ! Inorganic chemistry parameters
562      !----------------------------------
563      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909)
564909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' )
565
566      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910)
567910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' )
568
569      IF (lwp) THEN
570         WRITE(numsed,*) ' namelist nam_inorg'
571         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal
572         WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef
573         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal
574         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc
575         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc
576         WRITE(numsed,*) ' '
577      ENDIF
578
579      ! Unity conversion to get saturation conc. psat in [mol.l-1]
580      ! and reactivity rc in  [l.mol-1.s-1]
581      !----------------------------------------------------------
582      reac_sil   = rcopal / ryear     
583      reac_ligc  = rcligc / ryear
584
585      ! Additional parameter linked to POC/O2/No3/Po4
586      !----------------------------------------------
587      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911)
588911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' )
589
590      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912)
591912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' )
592
593      IF (lwp) THEN
594         WRITE(numsed,*) ' namelist nam_poc'
595         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2
596         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3
597         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4
598         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC
599         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep
600         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl
601         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs
602         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr
603         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4
604         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s
605         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2
606         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s
607         WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes
608         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso
609         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2
610         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3
611         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo
612         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4
613         WRITE(numsed,*) ' '
614      ENDIF
615
616
617      so2ut  = redO2    / redC
618      srno3  = redNo3   / redC
619      spo4r  = redPo4   / redC
620      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC
621      ! reactivity rc in  [l.mol-1.s-1]
622      reac_pocl  = rcorgl / ryear
623      reac_pocs  = rcorgs / ryear
624      reac_pocr  = rcorgr / ryear
625      reac_nh4   = rcnh4  / ryear
626      reac_h2s   = rch2s  / ryear
627      reac_fe2   = rcfe2  / ryear
628      reac_feh2s = rcfeh2s/ ryear
629      reac_fes   = rcfes  / ryear
630      reac_feso  = rcfeso / ryear
631
632      ! reactivity rc in  [l.mol-1.s-1]     
633      reac_cal = rccal / ryear
634
635      ! Bioturbation parameter
636      !------------------------
637      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913)
638913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' )
639
640      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914)
641914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' )
642
643      IF (lwp) THEN
644         WRITE(numsed,*) ' namelist nam_btb ' 
645         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot
646         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz
647         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc
648         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4
649         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig
650         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc
651         WRITE(numsed,*) ' '
652      ENDIF
653
654      ! Initial value (t=0) for sediment pore water and solid components
655      !----------------------------------------------------------------
656      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915)
657915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' )
658
659      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916)
660916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' )
661
662      IF (lwp) THEN
663         WRITE(numsed,*) ' namelist  nam_rst ' 
664         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
665         WRITE(numsed,*) ' '
666      ENDIF
667      nn_dtsed = 1
668
669
670   END SUBROUTINE sed_init_nam
671
672END MODULE sedini
Note: See TracBrowser for help on using the repository browser.