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/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedini.F90 @ 15297

Last change on this file since 15297 was 15297, checked in by aumont, 2 years ago

major update of the sediment module
  • Property svn:keywords set to Id
File size: 31.5 KB
Line 
1MODULE sedini
2   !!======================================================================
3   !!              ***  MODULE  sedini  ***
4   !! Sediment : define sediment variables
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   sed_ini    : 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), PUBLIC :: sedmask 
26
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
41   REAL(wp), PUBLIC    ::  &
42      redO2    =  138.  ,  &  !: Redfield coef for Oxygen
43      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate
44      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate
45      redC     =  117.  ,  &  !: 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      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1]
55      rcfesp   =   5E-6 ,  &  !: Precipitation of FeS [mol/l-1.an-1]
56      rcfesd   =   1E-3 ,  &  !: Dissolution of FeS [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/ 1.104E-5, 9.78E-6, 3.58E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6, 4.81E-6, 4.81E-6, 4.81E-6, 4.59E-6 /
74
75
76   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2
77   DATA diff2/ 4.47E-7, 3.89E-7, 1.77E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7, 2.51E-7, 2.51E-7, 2.51E-7, 1.74E-7 /
78
79   !! *  Routine accessibility
80   PUBLIC sed_ini          ! routine called by opa.F90
81
82   !! * Substitutions
83#  include "do_loop_substitute.h90"
84   !! $Id$
85CONTAINS
86
87
88   SUBROUTINE sed_ini
89      !!----------------------------------------------------------------------
90      !!                   ***  ROUTINE sed_ini  ***
91      !!
92      !! ** Purpose :  Initialization of sediment module
93      !!               - Reading namelist
94      !!               - Read the deepest water layer thickness
95      !!                 ( using as mask ) in Netcdf file
96      !!               - Convert unity if necessary
97      !!               - sets initial sediment composition
98      !!                 ( only clay or reading restart file )
99      !!               - sets sediment grid, porosity and others constants
100      !!
101      !!   History :
102      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
103      !!        !  06-07  (C. Ethe)  Re-organization
104      !!----------------------------------------------------------------------
105      INTEGER  :: ji, jj, js, jn, jk, ikt, ierr
106      REAL(wp) :: ztmp1, ztmp2 
107      !!----------------------------------------------------------------------
108
109      ! Reading namelist.sed variables
110      !---------------------------------------
111
112      CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
113
114      IF (lwp) THEN
115         WRITE(numsed,*)
116         WRITE(numsed,*) '                 PISCES framework'
117         WRITE(numsed,*) '                 SEDIMENT model'
118         WRITE(numsed,*) '                version 3.0  (2018) '
119         WRITE(numsed,*)
120         WRITE(numsed,*)
121      ENDIF
122
123      IF(lwp) WRITE(numsed,*) ' sed_ini : Initialization of sediment module  '
124      IF(lwp) WRITE(numsed,*) ' '
125
126      ! Read sediment Namelist
127      !-------------------------
128      CALL sed_ini_nam
129
130      ! Allocate SEDIMENT arrays
131      ierr =        sed_alloc()
132      ierr = ierr + sed_oce_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      gdepbot(:,:) = 0.
138      DO_2D( 0, 0, 0, 0 )
139         ikt = mbkt(ji,jj) 
140         IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ji,jj,ikt)
141         gdepbot(ji,jj) = gdepw_0(ji,jj,ikt+1)
142      END_2D
143
144      ! computation of total number of ocean points
145      !--------------------------------------------
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 )
153
154      ! Allocate memory size of global variables
155      ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) )
156      ALLOCATE( pwcpa(jpoce,jpksed,jpwat) )  ;  ALLOCATE( solcpa(jpoce,jpksed,jpsol) )
157      ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( rainrm_dta(jpoce,jpsol) )
158      ALLOCATE( rainrg(jpoce,jpsol) )
159      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) )
160      ALLOCATE( slatit(jpoce) )              ;  ALLOCATE( slongit(jpoce) )
161      ALLOCATE( zkbot(jpoce) )               ;  ALLOCATE( db(jpoce,jpksed) )
162      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) ) 
163      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) )
164      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) )
165      ALLOCATE( densSW(jpoce) )   ;   ALLOCATE( saturco3(jpoce,jpksed) ) 
166      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) )
167      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) )
168      ALLOCATE( rearatpom(jpoce, jpksed) )   ;  ALLOCATE( volc(jpoce,jpksed,jpsol) )
169      ALLOCATE( radssol(jpksed, jpwat) )     ;  ALLOCATE( rads1sol(jpksed, jpwat) )
170      ALLOCATE( apluss(jpoce, jpksed) )      ;  ALLOCATE( aminuss(jpoce,jpksed) )
171
172      ! Initialization of global variables
173      pwcp  (:,:,:) = 0.   ;  pwcp_dta  (:,:) = 0. 
174      pwcpa (:,:,:) = 0.   ;  solcpa(:,:,:) = 0.
175      solcp (:,:,:) = 0.   ;  rainrm_dta(:,:) = 0.
176      rainrg(:,:  ) = 0.
177      dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0.
178      temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0.
179      densSW (:   ) = 0.   ;  db     (:,:) = 0. 
180      hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0. 
181      dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0. 
182      fecratio(:)   = 1E-5 ;  rearatpom(:,:) = 0. 
183      radssol(:,:)  = 1.0  ;  rads1sol(:,:) = 0.
184      apluss(:,:)   = 0.0  ;  aminuss(:,:)  = 0.0
185
186      ! Chemical variables     
187      ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) )     
188      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )   
189      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )   
190      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) ) 
191      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) )
192      ALLOCATE( fluorids(jpoce) ) ; ALLOCATE( akh2s(jpoce) )    ;  ALLOCATE( aknh3(jpoce) )
193
194      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0.
195      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0.
196      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0.
197      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0. ;   akh2s  (:) = 0.
198      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0.
199      aknh3(:)  = 0.
200
201      ! Mass balance calculation 
202      ALLOCATE( fromsed(jpoce, jpsol+jpads) ) ; ALLOCATE( tosed(jpoce, jpsol+jpads) )
203
204      fromsed(:,:) = 0.    ;   tosed(:,:) = 0.
205
206      ! Initialization of sediment geometry
207      !------------------------------------
208      CALL sed_ini_geom
209
210      ! Offline specific mode
211      ! ---------------------
212      ln_sediment_offline = .FALSE.
213
214      !---------------------------------------------
215      ! Molecular weight [g/mol] for solid species
216      !---------------------------------------------
217
218      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
219      !---------------------------------------
220      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. )
221
222      !  clay
223      !  some kind of Illit (according to Pape)
224      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]
225      !--------------------------------------------------------------------
226      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
227         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    &
228         &              0.59 * 27. + 10. * 16.
229
230      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0)
231
232
233      mol_wgt(jsfes)  = 55.0 + 32.0
234
235      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
236      ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
237      ! So C(122)H(355)O(120)N(16)P(1)
238      !------------------------------------------------------------
239      mol_wgt(jspoc) = ( redC * 12. + 355. + 120. * 16. + redNo3 * 14. + 31. ) / redC
240      mol_wgt(jspos) = mol_wgt(jspoc)
241      mol_wgt(jspor) = mol_wgt(jspoc)
242
243      ! CaCO3
244      !---------
245      mol_wgt(jscal) = 40. + 12. + 3. * 16.
246
247      ! Density of solid material in sediment [g/cm**3]
248      !------------------------------------------------
249      ALLOCATE( dens_sol(jpsol) )
250      dens_sol(jsclay) = 2.6
251      dens_sol(jscal)  = 2.7
252      dens_sol(jsopal) = 2.1 
253      dens_sol(jsfeo)  = 3.4
254      dens_sol(jsfes)  = 4.8
255      dens_sol(jspoc)  = 1.0
256      dens_sol(jspos)  = 1.0
257      dens_sol(jspor)  = 1.0
258
259      ! Accumulation rate from Burwicz et al. (2011). This is used to
260      ! compute the flux of clays and minerals
261      ! --------------------------------------------------------------
262      DO ji = 1, jpoce
263          ztmp1 = 0.0117 * 3.0 / ( 1.0 + ( zkbot(ji) / 200.)**3 )
264          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 5000.)**10 )
265          wacc(ji) = ztmp2+ztmp1
266      END DO
267
268      ! Vertical profile of of the adsorption factor for adsorbed species
269      ! -----------------------------------------------------------------
270      radssol(:,jwfe2) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) )
271      radssol(:,jwnh4) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) )
272      rads1sol(:,jwnh4) = adsnh4 * por1(:) / ( por(:) + adsnh4 * por1(:) )
273      rads1sol(:,jwfe2) = adsfe2 * por1(:) / ( por(:) + adsfe2 * por1(:) )
274
275      ! Initialization of the array for non linear solving
276      ! --------------------------------------------------
277
278      ALLOCATE( jarr(jpksed*jpvode,2) )
279      ALLOCATE( jsvode(jpvode), isvode(jptrased) )
280      jsvode(1) = jwoxy ; jsvode(2) = jwno3 ; jsvode(3) = jwnh4
281      jsvode(4) = jwh2s ; jsvode(5) = jwso4 ; jsvode(6) = jwfe2
282      jsvode(7) = jpwat+jsfeo ; jsvode(8) = jpwat+jsfes
283      isvode(jwoxy) = 1 ; isvode(jwno3) = 2 ; isvode(jwnh4) = 3
284      isvode(jwh2s) = 4 ; isvode(jwso4) = 5 ; isvode(jwfe2) = 6
285      isvode(jpwat+jsfeo) = 7 ; isvode(jpwat+jsfes) = 8
286      DO js = 1, jpvode
287         DO jk = 1, jpksed
288            jn = (jk-1) * jpvode + js
289            jarr(jn,1) = jk
290            jarr(jn,2) = jsvode(js)
291         END DO
292      END DO
293
294      ALLOCATE( rstepros(jpoce) )
295
296#if defined key_sed_off
297      ln_sediment_offline = .TRUE.
298      IF (lwp) write(numsed,*) 'Sediment module is run in offline mode'
299      IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time'
300      IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false'
301      IF (lwp) write(numsed,*) '--------------------------------------------'
302      ln_sed_2way = .FALSE.
303#endif
304      ! Initialisation of tracer damping
305      ! --------------------------------
306      IF (ln_sediment_offline) THEN
307         CALL trc_dmp_sed_ini
308      ENDIF
309
310   END SUBROUTINE sed_ini
311
312   SUBROUTINE sed_ini_geom
313      !!----------------------------------------------------------------------
314      !!                   ***  ROUTINE sed_ini_geom  ***
315      !!
316      !! ** Purpose :  Initialization of sediment geometry
317      !!               - Read the deepest water layer thickness
318      !!                 ( using as mask ) in Netcdf file
319      !!               - sets sediment grid, porosity and molecular weight
320      !!                 and others constants
321      !!
322      !!   History :
323      !!        !  06-07  (C. Ethe)  Original
324      !!----------------------------------------------------------------------
325      !! * Modules used
326      !! * local declarations
327      INTEGER  :: ji, jj, jk, jn
328      REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw
329      REAL(wp) :: ztmp1, ztmp2
330      !----------------------------------------------------------
331
332      IF(lwp) WRITE(numsed,*) ' sed_ini_geom : Initialization of sediment geometry '
333      IF(lwp) WRITE(numsed,*) ' '
334
335      ! Computation of 1D array of sediments points
336      indoce = 0
337      DO_2D( 0, 0, 0, 0 )
338         IF (  epkbot(ji,jj) > 0. ) THEN
339            indoce          = indoce + 1
340            iarroce(indoce) = (jj - 1) * jpi + ji
341         ENDIF
342      END_2D
343
344      IF ( indoce .EQ. 0 ) THEN
345         indoce = 1
346         iarroce(indoce) = 1
347      ENDIF
348
349      IF( indoce .NE. jpoce ) THEN
350         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' )
351      ELSE
352         IF (lwp) WRITE(numsed,*) ' '
353         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
354         IF (lwp) WRITE(numsed,*) ' '
355      ENDIF
356
357      ! initialization of dzkbot in [cm]
358      !------------------------------------------------   
359      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
360#if defined key_sed_off
361      dzkbot(1:jpoce) = 1.E8
362#else
363      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
364#endif
365      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) )
366      CALL pack_arr ( jpoce, slatit(1:jpoce), gphit(1:jpi,1:jpj), iarroce(1:jpoce) )
367      CALL pack_arr ( jpoce, slongit(1:jpoce), glamt(1:jpi,1:jpj), iarroce(1:jpoce) )
368
369      ! Geometry and  constants
370      ! sediment layer thickness [cm]
371      ! (1st layer= diffusive layer = pur water)
372      !------------------------------------------
373      za1  = (  sedzmin - sedhmax / FLOAT(jpksed-1)  )                                                      &
374         & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * (  LOG( COSH( (jpksed - sedkth) / sedacr) )      &
375         &                                                   - LOG( COSH( ( 1  - sedkth) / sedacr) )  )  )
376      za0  = sedzmin - za1 * TANH( (1-sedkth) / sedacr )
377      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  )
378
379      dz(1)       = 0.1
380      profsedw(1) = 0.0
381      profsed(1)  = -dz(1) / 2.
382      DO jk = 2, jpksed
383         zw = REAL( jk , wp )
384         zt = REAL( jk , wp ) - 0.5_wp
385         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  ) 
386         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  )
387         dz(jk) = profsedw(jk) - profsedw(jk-1)
388      END DO
389
390      DO ji = 1, jpoce
391         dz3d(ji,:) = dz(:)
392      END DO
393
394      !  Porosity profile [0]
395      !---------------------
396      por(1) = 1.0
397      DO jk = 2, jpksed
398         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * profsed(jk) )
399      END DO
400 
401      ! inverse of  Porosity profile
402      !-----------------------------
403      por1(:) = 1. - por(:)
404
405      ! Volumes of pore water and solid fractions (vector and array)
406      !     WARNING : volw(1) and vols(1) are sublayer volums
407      volw(:) = dz(:) * por(:)
408      vols(:) = dz(:) * por1(:)
409
410      ! temporary new value for dz3d(:,1)
411      dz3d(1:jpoce,1) = dzkbot(1:jpoce)
412
413      ! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
414      DO jk = 1, jpksed
415         volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
416         vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
417      ENDDO
418
419      ! Back to the old sublayer vlaue for dz3d(:,1)
420      dz3d(1:jpoce,1) = dz(1)
421
422   END SUBROUTINE sed_ini_geom
423
424   SUBROUTINE sed_ini_nam
425      !!----------------------------------------------------------------------
426      !!                   ***  ROUTINE sed_ini_nam  ***
427      !!
428      !! ** Purpose :  Initialization of sediment geometry
429      !!               - Reading namelist and defines constants variables
430      !!
431      !!   History :
432      !!        !  06-07  (C. Ethe)  Original
433      !!----------------------------------------------------------------------
434
435      INTEGER ::   numnamsed_ref = -1           !! Logical units for namelist sediment
436      INTEGER ::   numnamsed_cfg = -1           !! Logical units for namelist sediment
437      INTEGER :: ios                 ! Local integer output status for namelist read
438      CHARACTER(LEN=20)   ::   clname
439
440      TYPE PSED
441         CHARACTER(len = 20)  :: snamesed   !: short name
442         CHARACTER(len = 80 ) :: lnamesed   !: long name
443         CHARACTER(len = 20 ) :: unitsed    !: unit
444      END TYPE PSED
445
446      TYPE(PSED) , DIMENSION(jpsol     ) :: sedsol
447      TYPE(PSED) , DIMENSION(jpwat     ) :: sedwat
448      TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
449      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
450
451      NAMELIST/nam_run/ln_sed_2way,nrosorder,rosatol,rosrtol
452      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox
453      NAMELIST/nam_trased/sedsol, sedwat
454      NAMELIST/nam_diased/seddiag3d, seddiag2d
455      NAMELIST/nam_inorg/rcopal, rccal, ratligc, rcligc
456      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  &
457         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfeso, rcfesp, &
458         &             rcfesd, xksedo2, xksedno3, xksedfeo, xksedso4
459      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, adsfe2, ln_irrig, xirrzsc
460      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out
461
462      INTEGER :: ji, jn, jn1
463      !-------------------------------------------------------
464
465      IF(lwp) WRITE(numsed,*) ' sed_ini_nam : Read namelists '
466      IF(lwp) WRITE(numsed,*) ' '
467
468      ! ryear = 1 year converted in second
469      !------------------------------------
470      IF (lwp) THEN
471         WRITE(numsed,*) ' '
472         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
473         WRITE(numsed,*) ' '     
474      ENDIF
475
476      ! Reading namelist.sed variables
477      !---------------------------------
478      clname = 'namelist_sediment'
479      IF(lwp) WRITE(numsed,*) ' sed_ini_nam : read SEDIMENT namelist'
480      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~'
481      CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
482      CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
483
484      nitsed000 = nittrc000
485      nitsedend = nitend
486      ! Namelist nam_run
487      REWIND( numnamsed_ref )              ! Namelist nam_run in reference namelist : Pisces variables
488      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901)
489901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' )
490
491      REWIND( numnamsed_cfg )              ! Namelist nam_run in reference namelist : Pisces variables
492      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902)
493902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' )
494
495      IF (lwp) THEN
496         WRITE(numsed,*) ' namelist nam_run'
497         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way
498         WRITE(numsed,*) ' Order of the Rosenbrock method (2,3,4) = ', nrosorder
499         WRITE(numsed,*) ' Tolerance for absolute error = ', rosatol
500         WRITE(numsed,*) ' Tolerance for relative order = ', rosrtol
501      ENDIF
502
503      IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' )
504
505      REWIND( numnamsed_ref )              ! Namelist nam_geom in reference namelist : Pisces variables
506      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
507903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' )
508
509      REWIND( numnamsed_cfg )              ! Namelist nam_geom in reference namelist : Pisces variables
510      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904)
511904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' )
512
513      IF (lwp) THEN
514         WRITE(numsed,*) ' namelist nam_geom'
515         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed
516         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin
517         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax
518         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth
519         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr
520         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf
521         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf
522         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox
523      ENDIF
524
525      jpksedm1  = jpksed - 1
526      dtsed = rDt_trc
527
528      REWIND( numnamsed_ref )              ! Namelist nam_trased in reference namelist : Pisces variables
529      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905)
530905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' )
531
532      REWIND( numnamsed_cfg )              ! Namelist nam_trased in reference namelist : Pisces variables
533      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906)
534906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' )
535
536      DO jn = 1, jpsol
537         sedtrcd(jn) = sedsol(jn)%snamesed
538         sedtrcl(jn) = sedsol(jn)%lnamesed
539         sedtrcu(jn) = sedsol(jn)%unitsed
540      END DO
541
542      DO jn = 1, jpwat
543         jn1 = jn + jpsol
544         sedtrcd(jn1) = sedwat(jn)%snamesed
545         sedtrcl(jn1) = sedwat(jn)%lnamesed
546         sedtrcu(jn1) = sedwat(jn)%unitsed
547      END DO
548
549      IF (lwp) THEN
550         WRITE(numsed,*) ' namelist nam_trased'
551         WRITE(numsed,*) ' '
552         DO jn = 1, jptrased
553            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
554            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
555            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
556            WRITE(numsed,*) ' '
557         END DO
558         WRITE(numsed,*) ' '
559      ENDIF
560
561      REWIND( numnamsed_ref )              ! Namelist nam_diased in reference namelist : Pisces variables
562      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907)
563907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' )
564
565      REWIND( numnamsed_cfg )              ! Namelist nam_diased in reference namelist : Pisces variables
566      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908)
567908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' )
568     
569      DO jn = 1, jpdia3dsed
570         seddia3d(jn) = seddiag3d(jn)%snamesed
571         seddia3l(jn) = seddiag3d(jn)%lnamesed
572         seddia3u(jn) = seddiag3d(jn)%unitsed
573      END DO
574
575      DO jn = 1, jpdia2dsed
576         seddia2d(jn) = seddiag2d(jn)%snamesed
577         seddia2l(jn) = seddiag2d(jn)%lnamesed
578         seddia2u(jn) = seddiag2d(jn)%unitsed
579      END DO
580
581      IF (lwp) THEN
582         WRITE(numsed,*) ' namelist nam_diased'
583         WRITE(numsed,*) ' '
584         DO jn = 1, jpdia3dsed
585            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
586            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
587            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
588            WRITE(numsed,*) ' '
589         END DO
590
591         DO jn = 1, jpdia2dsed
592            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
593            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
594            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
595            WRITE(numsed,*) ' '
596         END DO
597
598         WRITE(numsed,*) ' '
599      ENDIF
600
601      ! Inorganic chemistry parameters
602      !----------------------------------
603      REWIND( numnamsed_ref )              ! Namelist nam_inorg in reference namelist : Pisces variables
604      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909)
605909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' )
606
607      REWIND( numnamsed_cfg )              ! Namelist nam_inorg in reference namelist : Pisces variables
608      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910)
609910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' )
610
611      IF (lwp) THEN
612         WRITE(numsed,*) ' namelist nam_inorg'
613         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal
614         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal
615         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc
616         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc
617         WRITE(numsed,*) ' '
618      ENDIF
619
620      ! Unity conversion to get saturation conc. psat in [mol.l-1]
621      ! and reactivity rc in  [l.mol-1.s-1]
622      !----------------------------------------------------------
623      reac_sil   = rcopal / ryear     
624      reac_ligc  = rcligc / ryear
625
626      ! Additional parameter linked to POC/O2/No3/Po4
627      !----------------------------------------------
628      REWIND( numnamsed_ref )              ! Namelist nam_poc in reference namelist : Pisces variables
629      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911)
630911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' )
631
632      REWIND( numnamsed_cfg )              ! Namelist nam_poc in reference namelist : Pisces variables
633      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912)
634912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' )
635
636      IF (lwp) THEN
637         WRITE(numsed,*) ' namelist nam_poc'
638         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2
639         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3
640         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4
641         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC
642         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep
643         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl
644         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs
645         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr
646         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4
647         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s
648         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2
649         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s
650         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso
651         WRITE(numsed,*) ' Precipitation of FeS             rcfesp   = ', rcfesp
652         WRITE(numsed,*) ' Dissolution of FeS               rcfesd   = ', rcfesd
653         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2
654         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3
655         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo
656         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4
657         WRITE(numsed,*) ' '
658      ENDIF
659
660
661      so2ut  = redO2    / redC
662      srno3  = redNo3   / redC
663      spo4r  = redPo4   / redC
664      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC
665      ! reactivity rc in  [l.mol-1.s-1]
666      reac_pocl  = rcorgl / ryear
667      reac_pocs  = rcorgs / ryear
668      reac_pocr  = rcorgr / ryear
669      reac_nh4   = rcnh4  / ryear
670      reac_h2s   = rch2s  / ryear
671      reac_fe2   = rcfe2  / ryear
672      reac_feh2s = rcfeh2s/ ryear
673      reac_feso  = rcfeso / ryear
674      reac_fesp  = rcfesp / ryear
675      reac_fesd  = rcfesd / ryear
676
677      ! reactivity rc in  [l.mol-1.s-1]     
678      reac_cal = rccal / ryear
679
680      ! Bioturbation parameter
681      !------------------------
682      REWIND( numnamsed_ref )              ! Namelist nam_btb in reference namelist : Pisces variables
683      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913)
684913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' )
685
686      REWIND( numnamsed_cfg )              ! Namelist nam_btb in reference namelist : Pisces variables
687      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914)
688914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' )
689
690      IF (lwp) THEN
691         WRITE(numsed,*) ' namelist nam_btb ' 
692         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot
693         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz
694         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc
695         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4
696         WRITE(numsed,*) ' Adsorption coefficient of Fe2     adsfe2   = ', adsfe2
697         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig
698         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc
699         WRITE(numsed,*) ' '
700      ENDIF
701
702      ! Initial value (t=0) for sediment pore water and solid components
703      !----------------------------------------------------------------
704      REWIND( numnamsed_ref )              ! Namelist nam_rst in reference namelist : Pisces variables
705      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915)
706915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' )
707
708      REWIND( numnamsed_cfg )              ! Namelist nam_rst in reference namelist : Pisces variables
709      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916)
710916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' )
711
712      IF (lwp) THEN
713         WRITE(numsed,*) ' namelist  nam_rst ' 
714         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
715         WRITE(numsed,*) ' '
716      ENDIF
717
718      CLOSE( numnamsed_cfg )
719      CLOSE( numnamsed_ref )
720
721   END SUBROUTINE sed_ini_nam
722
723END MODULE sedini
Note: See TracBrowser for help on using the repository browser.