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 trunk/NEMO/TOP_SRC/SED – NEMO

source: trunk/NEMO/TOP_SRC/SED/sedini.F90 @ 1264

Last change on this file since 1264 was 1264, checked in by cetlod, 15 years ago

clean TOP model routines to avoid warning when compiling, see ticket:303

File size: 34.2 KB
Line 
1MODULE sedini
2   !!======================================================================
3   !!              ***  MODULE  sedini  ***
4   !! Sediment : define sediment variables
5   !!=====================================================================
6#if defined key_sed
7
8   !!----------------------------------------------------------------------
9   !!   sed_init    : initialization, namelist read, and parameters control
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE sed     ! sediment global variable
13   USE seddta
14   USE sedrst
15   USE sedco3
16   USE sedchem
17   USE sedarr
18   USE iom
19
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! Module variables
25   REAL(wp)    ::  &
26      sisat   =  800.  ,  &  !: saturation for si    [ mol.l-1]
27      claysat =    0.  ,  &  !: saturation for clay  [ mol.l-1]
28      rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1]
29      rcclay  =    0.  ,  &  !: reactivity for clay  [l.mol-1.an-1]
30      dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s]
31
32   REAL(wp)    ::  &
33      redO2   =  172.  ,  &  !: Redfield coef for Oxygene
34      redNo3  =   16.  ,  &  !: Redfield coef for Nitrates
35      redPo4  =    1.  ,  &  !: Redfield coef for Phosphates
36      redC    =  122.  ,  &  !: Redfield coef for Carbone
37      redDnit =  97.6  ,  &  !: Redfield coef for denitrification   
38      rcorg   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]   
39      o2seuil =    1.  ,  &  !: threshold O2 concentration for [mol.l-1]     
40      rcorgN  =   50.       !: reactivity for POC/No3 [l.mol-1.an-1]
41
42   REAL(wp)    ::  &
43      rccal   = 1000.        !: reactivity for calcite         [l.mol-1.an-1]
44
45   REAL(wp)    ::  &
46      dbiot   = 15.          !: coefficient for bioturbation    [cm**2.(1000an-1)]
47
48   LOGICAL     :: &
49      ln_rst_sed = .TRUE.    !: initialisation from a restart file or not
50
51   REAL(wp)    ::  &
52      ryear = 365. * 24. * 3600. !:  1 year converted in second
53
54   !! *  Routine accessibility
55   PUBLIC sed_init          ! routine called by opa.F90
56
57CONTAINS
58
59
60   SUBROUTINE sed_init
61      !!----------------------------------------------------------------------
62      !!                   ***  ROUTINE sed_init  ***
63      !!
64      !! ** Purpose :  Initialization of sediment module
65      !!               - Reading namelist
66      !!               - Read the deepest water layer thickness
67      !!                 ( using as mask ) in Netcdf file
68      !!               - Convert unity if necessary
69      !!               - sets initial sediment composition
70      !!                 ( only clay or reading restart file )
71      !!               - sets sediment grid, porosity and others constants
72      !!
73      !!   History :
74      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
75      !!        !  06-07  (C. Ethe)  Re-organization
76      !!----------------------------------------------------------------------
77      INTEGER :: ji, jj, ikt
78      CHARACTER (len=28) :: csedname 
79#if defined key_sed_off
80      INTEGER  :: numblt         
81      INTEGER  :: nummsh   
82      REAL(wp) , DIMENSION(jpi,jpj) :: zdta
83#endif
84      !!----------------------------------------------------------------------
85
86
87      ! Reading namelist.sed variables
88      !---------------------------------------
89
90      csedname = 'sediment.output'
91      CALL ctlopn( numsed, csedname, 'OLD', 'FORMATTED', 'SEQUENTIAL', 1, numout, .FALSE., 1 )
92
93      WRITE(numsed,*)
94      WRITE(numsed,*) '                 L S C E - C E A'
95      WRITE(numsed,*) '                 SEDIMENT model'
96      WRITE(numsed,*) '                version 2.0  (2007) '
97      WRITE(numsed,*)
98      WRITE(numsed,*)
99
100   
101      WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
102      WRITE(numsed,*) ' '
103
104      ! Determination of sediments number of points and allocate global variables
105#if defined key_sed_off
106      ! Reading Netcdf Pisces file for deepest water layer thickness [m]
107      ! This data will be used as mask to merdge space in 1D vector
108      !----------------------------------------------------------------
109
110      CALL iom_open ( 'mesh_mask'  , nummsh )   
111      CALL iom_open ( 'e3tbot'     , numblt )
112
113      ! mask sediment points for outputs
114      CALL iom_get( nummsh, jpdom_data, 'tmask' , tmask )
115      CALL iom_get( nummsh, jpdom_data, 'mbathy', sbathy )
116     
117      ! longitude/latitude values
118      CALL iom_get ( nummsh, jpdom_data,'nav_lon', glamt(:,:) )
119      CALL iom_get ( nummsh, jpdom_data,'nav_lat', gphit(:,:) )
120   
121      ! bottom layer thickness
122      CALL iom_get  ( numblt, jpdom_data, 'E3TBOT', zdta(:,:) )
123      epkbot(:,:) = 0.
124      DO jj = 1, jpj
125         DO ji = 1, jpi
126            ikt = MAX( INT( sbathy(ji,jj) )  - 1, 1 )
127            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = zdta(ji,jj)
128         ENDDO
129      ENDDO
130
131      CALL iom_close( nummsh ) 
132      CALL iom_close( numblt ) 
133#else
134
135      epkbot(:,:) = 0.
136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            ikt = MAX( mbathy(ji,jj) - 1, 1 )
139            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ikt)
140         ENDDO
141      ENDDO
142#endif     
143
144
145      ! computation of total number of ocean points
146      !--------------------------------------------
147      jpoce  = COUNT( epkbot(:,:) > 0. )
148
149
150      ! Allocate memory size of global variables
151      ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) )
152      ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ;  ALLOCATE( rainrm_dta(jpoce,jpsol) )
153      ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) ) 
154      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) )
155      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )               
156      ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) ) 
157      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) ) 
158      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) )
159
160      ! Initialization of global variables
161      pwcp  (:,:,:) = 0.  ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0. 
162      solcp (:,:,:) = 0.  ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0.
163      rainrm(:,:  ) = 0.  ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0. 
164      dzdep (:    ) = 0.  ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0.
165      temp  (:    ) = 0.  ;  salt   (:   ) = 0. 
166      press (:    ) = 0.  ;  densSW (:   ) = 0. 
167      hipor (:,:  ) = 0.  ;  co3por (:,: ) = 0. 
168      dz3d  (:,:  ) = 0.  ;  volw3d (:,: ) = 0.  ;  vols3d   (:,:) = 0. 
169
170      ! Chemical variables     
171      ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) )     
172      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )   
173      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )   
174      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )   
175
176      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0.
177      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0.
178      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0.
179      borats(:) = 0. ;   calcon2(:) = 0.
180
181      ! Mass balance calculation 
182      ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ;  ALLOCATE( rloss(jpoce, jpsol) )
183      ALLOCATE( tokbot (jpoce, jpwat) ) 
184
185      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0. 
186
187      ! Read sediment Namelist
188      !-------------------------
189      CALL sed_init_nam
190   
191      ! Initialization of sediment geometry
192      !------------------------------------
193      CALL sed_init_geom
194
195
196      ! sets initial sediment composition
197      ! ( only clay or reading restart file )
198      !---------------------------------------
199      CALL sed_init_data
200
201
202      CALL sed_init_wri
203
204
205   END SUBROUTINE sed_init
206
207
208   SUBROUTINE sed_init_geom
209      !!----------------------------------------------------------------------
210      !!                   ***  ROUTINE sed_init_geom  ***
211      !!
212      !! ** Purpose :  Initialization of sediment geometry
213      !!               - Read the deepest water layer thickness
214      !!                 ( using as mask ) in Netcdf file
215      !!               - sets sediment grid, porosity and molecular weight
216      !!                 and others constants
217      !!
218      !!   History :
219      !!        !  06-07  (C. Ethe)  Original
220      !!----------------------------------------------------------------------
221      !! * Modules used
222      !! * local declarations
223      INTEGER ::  &
224        ji, jj, jk
225 
226#if defined key_sed_off
227      REAL(wp) , DIMENSION (jpi,jpj) :: zdta
228      INTEGER  ::  numpres = 52
229#endif
230      !----------------------------------------------------------
231
232      WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
233      WRITE(numsed,*) ' '
234
235      ! Computation of 1D array of sediments points
236      indoce = 0
237      DO jj = 1, jpj
238         DO ji = 1, jpi
239            IF (  epkbot(ji,jj) > 0. ) THEN
240               indoce          = indoce + 1
241               iarroce(indoce) = (jj - 1) * jpi + ji
242            ENDIF
243         END DO
244      END DO
245
246      IF( indoce .NE. jpoce ) THEN
247         WRITE(numsed,*) ' '
248         WRITE(numsed,*) 'Warning : number of ocean points indoce = ', indoce, &
249            &     ' doesn''t match  number of point where epkbot>0  jpoce = ', jpoce
250         WRITE(numsed,*) ' '
251         WRITE(numsed,*) ' '
252         STOP
253      ELSE
254         WRITE(numsed,*) ' '
255         WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
256         WRITE(numsed,*) ' '
257      ENDIF
258
259#if defined key_sed_off
260
261      ! Reading Netcdf Pisces file for deepest water layer thickness [m]
262      ! This data will be used as mask to merdge space in 1D vector
263      !----------------------------------------------------------------
264      CALL iom_open ( 'pressbot'   , numpres )
265     
266      ! pressure  in  bars     
267      CALL iom_get ( numpres, jpdom_data,'BATH', zdta(:,:) )
268      CALL pack_arr ( jpoce, press(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) )
269      press(1:jpoce) = press(1:jpoce) * 1.025e-1 
270
271      CALL iom_close ( numpres )
272#endif
273
274
275      ! mask sediment points for outputs
276      DO jk = 1, jpksed
277         tmasksed(:,:,jk) = tmask(:,:,1)
278      ENDDO
279
280      ! initialization of dzkbot in [cm]
281      !------------------------------------------------   
282      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
283      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
284
285      ! Geometry and  constants
286      ! sediment layer thickness [cm]
287      ! (1st layer= diffusive layer = pur water)
288      !------------------------------------------
289      dz(1)  = 0.1
290      dz(2)  = 0.3
291      dz(3)  = 0.3
292      dz(4)  = 0.5
293      dz(5)  = 0.5
294      dz(6)  = 0.5
295      dz(7)  = 1.
296      dz(8)  = 1.
297      dz(9)  = 1.
298      dz(10) = 2.45
299      dz(11) = 2.45
300
301      DO jk = 1, jpksed
302         DO ji = 1, jpoce
303            dz3d(ji,jk) = dz(jk)
304         END DO
305      ENDDO
306
307      ! Depth(jk)= depth of middle of each layer
308      !----------------------------------------
309      profsed(1) = -dz(1)/ 2.
310      DO jk = 2, jpksed
311         profsed(jk) = profsed(jk-1) + dz(jk-1) / 2. + dz(jk) / 2.
312      ENDDO
313
314      !  Porosity profile [0]
315      !---------------------
316      por(1)  = 1.
317      por(2)  = 0.95
318      por(3)  = 0.9
319      por(4)  = 0.85
320      por(5)  = 0.81
321      por(6)  = 0.75
322      por(7)  = 0.75     
323      por(8)  = 0.75     
324      por(9)  = 0.75     
325      por(10) = 0.75
326      por(11) = 0.75
327 
328      ! inverse of  Porosity profile
329      !-----------------------------
330      por1(:) = 1. - por(:)
331
332      ! Volumes of pore water and solid fractions (vector and array)
333      !     WARNING : volw(1) and vols(1) are sublayer volums
334      volw(:) = dz(:) * por(:)
335      vols(:) = dz(:) * por1(:)
336
337      ! temporary new value for dz3d(:,1)
338      dz3d(1:jpoce,1) = dzkbot(1:jpoce)
339
340      ! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
341      DO jk = 1, jpksed
342         volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
343         vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
344      ENDDO
345
346      ! Back to the old sublayer vlaue for dz3d(:,1)
347      dz3d(1:jpoce,1) = dz(1)
348
349
350      !---------------------------------------------
351      ! Molecular weight [g/mol] for solid species
352      !---------------------------------------------
353
354
355      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
356      !---------------------------------------
357      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 
358
359      !  clay
360      !  some kind of Illit (according to Pape)
361      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]
362      !--------------------------------------------------------------------
363      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
364         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    &
365         &              0.59 * 27. + 10. * 16.
366
367      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
368      ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
369      ! So C(122)H(355)O(120)N(16)P(1)
370      !------------------------------------------------------------
371      mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
372         &                16. * 14. + 31. ) / 122.
373
374      ! CaCO3
375      !---------
376      mol_wgt(jscal) = 40. + 12. + 3. * 16.
377
378      ! Density of solid material in sediment [g/cm**3]
379      !------------------------------------------------
380      dens = 2.6
381
382      ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
383      !--------------------------------------------------------------------
384      diff(:) = dcoef * por(:)
385
386
387      ! Initialization of time step as function of porosity [cm**2/s]
388      !------------------------------------------------------------------
389      rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) )
390     
391   END SUBROUTINE sed_init_geom
392
393   SUBROUTINE sed_init_nam
394      !!----------------------------------------------------------------------
395      !!                   ***  ROUTINE sed_init_nam  ***
396      !!
397      !! ** Purpose :  Initialization of sediment geometry
398      !!               - Reading namelist and defines constants variables
399      !!
400      !!   History :
401      !!        !  06-07  (C. Ethe)  Original
402      !!----------------------------------------------------------------------
403
404      INTEGER :: &
405         numnamsed = 28
406
407      TYPE PSED
408         CHARACTER(len = 20)  :: snamesed   !: short name
409         CHARACTER(len = 80 ) :: lnamesed   !: long name
410         CHARACTER(len = 20 ) :: unitsed    !: unit
411      END TYPE PSED
412
413      TYPE(PSED) , DIMENSION(jpsol     ) :: sedsol
414      TYPE(PSED) , DIMENSION(jpwat     ) :: sedwat
415      TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
416      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
417
418      NAMELIST/nam_time/nfreq
419      NAMELIST/nam_trased/sedsol, sedwat
420      NAMELIST/nam_diased/seddiag3d, seddiag2d
421      NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef 
422      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, &
423         &           rcorg, o2seuil, rcorgN
424      NAMELIST/nam_cal/rccal
425      NAMELIST/nam_dc13/pdb, rc13P, rc13Ca
426      NAMELIST/nam_btb/dbiot
427      NAMELIST/nam_rst/ln_rst_sed
428
429      INTEGER :: jn, jn1
430
431      !-------------------------------------------------------
432
433      WRITE(numsed,*) ' sed_init_nam : Read namelists '
434      WRITE(numsed,*) ' '
435
436      ! ryear = 1 year converted in second
437      !------------------------------------
438      WRITE(numsed,*) ' '
439      WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
440      WRITE(numsed,*) ' '     
441
442      ! Reading namelist.sed variables
443      !---------------------------------------
444      OPEN( UNIT = numnamsed, FILE = 'namelist.sediment', FORM = 'FORMATTED', STATUS = 'OLD' )
445
446      dtsed = rdt
447#if ! defined key_sed_off
448      nitsed000 = nittrc000
449      nitsedend = nitend
450      nwrised   = nwritetrc
451#else
452      nitsed000 = nit000
453      nitsedend = nitend
454      nwrised   = nwrite
455#endif
456     ! Diffraction/reaction parameters
457      !----------------------------------
458      READ( numnamsed, nam_time )
459      WRITE(numsed,*) ' namelist nam_time'
460
461#if ! defined key_sed_off
462      nfreq     = 1
463#endif
464
465      WRITE(numsed,*) ' sedimentation time step              dtsed      = ', dtsed
466      WRITE(numsed,*) ' 1st time step for sediment.          nitsed000  = ', nitsed000
467      WRITE(numsed,*) ' last time step for sediment.         nitsedend  = ', nitsedend
468      WRITE(numsed,*) ' frequency of sediment outputs        nwrised    = ', nwrised
469      WRITE(numsed,*) ' frequency of restoring inputs data   nfreq      = ', nfreq
470      WRITE(numsed,*) ' '
471
472      REWIND( numnamsed )               ! read nattrc
473      READ  ( numnamsed, nam_trased )
474
475      DO jn = 1, jpsol
476         sedtrcd(jn) = sedsol(jn)%snamesed
477         sedtrcl(jn) = sedsol(jn)%lnamesed
478         sedtrcu(jn) = sedsol(jn)%unitsed
479      END DO
480
481      DO jn = 1, jpwat
482         jn1 = jn + jpsol
483         sedtrcd(jn1) = sedwat(jn)%snamesed
484         sedtrcl(jn1) = sedwat(jn)%lnamesed
485         sedtrcu(jn1) = sedwat(jn)%unitsed
486      END DO
487
488      WRITE(numsed,*) ' namelist nam_trased'
489      WRITE(numsed,*) ' '
490      DO jn = 1, jptrased
491         WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
492         WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
493         WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
494         WRITE(numsed,*) ' '
495      END DO
496      WRITE(numsed,*) ' '
497
498     
499      REWIND( numnamsed )
500      READ( numnamsed, nam_diased )
501
502      DO jn = 1, jpdia3dsed
503         seddia3d(jn) = seddiag3d(jn)%snamesed
504         seddia3l(jn) = seddiag3d(jn)%lnamesed
505         seddia3u(jn) = seddiag3d(jn)%unitsed
506      END DO
507
508      DO jn = 1, jpdia2dsed
509         seddia2d(jn) = seddiag2d(jn)%snamesed
510         seddia2l(jn) = seddiag2d(jn)%lnamesed
511         seddia2u(jn) = seddiag2d(jn)%unitsed
512      END DO
513
514      WRITE(numsed,*) ' namelist nam_diased'
515      WRITE(numsed,*) ' '
516      DO jn = 1, jpdia3dsed
517         WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
518         WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
519         WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
520         WRITE(numsed,*) ' '
521      END DO
522
523      DO jn = 1, jpdia2dsed
524         WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
525         WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
526         WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
527         WRITE(numsed,*) ' '
528      END DO
529
530      WRITE(numsed,*) ' '
531
532
533      ! Diffraction/reaction parameters
534      !----------------------------------
535      REWIND( numnamsed )
536      READ( numnamsed, nam_reac )
537      WRITE(numsed,*) ' namelist nam_reac'
538      WRITE(numsed,*) ' saturation for si    sisat   = ', sisat
539      WRITE(numsed,*) ' saturation for clay  claysat = ', claysat
540      WRITE(numsed,*) ' reactivity for Si    rcopal  = ', rcopal
541      WRITE(numsed,*) ' reactivity for clay  rcclay  = ', rcclay
542      WRITE(numsed,*) ' diff. coef for por.  dcoef   = ', dcoef
543      WRITE(numsed,*) ' '
544
545
546      ! Unity conversion to get saturation conc. psat in [mol.l-1]
547      ! and reactivity rc in  [l.mol-1.s-1]
548      !----------------------------------------------------------
549      sat_sil   = sisat   * 1.e-6   
550      sat_clay  = claysat * 1.e-6
551
552      reac_sil  = rcopal / ryear     
553      reac_clay  = rcclay / ryear
554
555
556      ! Additional parameter linked to POC/O2/No3/Po4
557      !----------------------------------------------
558      REWIND( numnamsed )
559      READ( numnamsed, nam_poc )
560      WRITE(numsed,*) ' namelist nam_poc'
561      WRITE(numsed,*) ' Redfield coef for oxy     redO2   = ', redO2
562      WRITE(numsed,*) ' Redfield coef for no3     redNo3  = ', redNo3
563      WRITE(numsed,*) ' Redfield coef for po4     redPo4  = ', redPo4
564      WRITE(numsed,*) ' Redfield coef for carbone redC    = ', redC
565      WRITE(numsed,*) ' Redfield coef for denitri redDnit = ', redDnit
566      WRITE(numsed,*) ' reactivity for POC/O2     rcorg   = ', rcorg
567      WRITE(numsed,*) ' threshold O2 concen.      o2seuil = ', o2seuil
568      WRITE(numsed,*) ' reactivity for POC/NO3    rcorgN  = ', rcorgN
569      WRITE(numsed,*) ' '
570
571
572      so2ut  = redO2   / redC
573      srno3  = redNo3  / redC
574      spo4r  = redPo4  / redC
575      srDnit = redDnit / redC
576      sthro2 = o2seuil * 1.e-6 ! threshold O2 concen. in [mol.l-1]
577      ! reactivity rc in  [l.mol-1.s-1]
578      reac_poc  = rcorg / ryear
579      reac_no3 = rcorgN  / ryear
580
581
582      ! Carbonate parameters
583      !---------------------
584      READ( numnamsed, nam_cal )
585      WRITE(numsed,*) ' namelist  nam_cal'
586      WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal
587      WRITE(numsed,*) ' '
588
589      ! reactivity rc in  [l.mol-1.s-1]     
590      reac_cal = rccal / ryear
591
592
593      ! C13  parameters
594      !----------------
595      READ( numnamsed, nam_dc13 )
596      WRITE(numsed,*) ' namelist nam_dc13 ' 
597      WRITE(numsed,*) ' 13C/12C in PD Belemnite        PDB    = ', pdb
598      WRITE(numsed,*) ' 13C/12C in POC = rc13P*PDB     rc13P  = ', rc13P
599      WRITE(numsed,*) ' 13C/12C in CaCO3 = rc13Ca*PDB  rc13Ca = ', rc13Ca
600      WRITE(numsed,*) ' '
601
602     
603      ! Bioturbation parameter
604      !------------------------
605      READ( numnamsed, nam_btb )
606      WRITE(numsed,*) ' namelist nam_btb ' 
607      WRITE(numsed,*) ' coefficient for bioturbation   dbiot    = ', dbiot
608      WRITE(numsed,*) ' '
609
610      ! Unity convertion to get bioturb coefficient in [cm2.s-1]
611      db = dbiot / ( ryear * 1000. )
612
613      ! Initial value (t=0) for sediment pore water and solid components
614      !----------------------------------------------------------------
615      READ( numnamsed, nam_rst )
616      WRITE(numsed,*) ' namelist  nam_rst ' 
617      WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
618      WRITE(numsed,*) ' '
619
620      CLOSE( numnamsed )
621
622   END SUBROUTINE sed_init_nam
623
624   SUBROUTINE sed_init_data
625      !!----------------------------------------------------------------------
626      !!                   ***  ROUTINE sed_init_data  ***
627      !!
628      !! ** Purpose :  Initialization of sediment module
629      !!               - sets initial sediment composition
630      !!                 ( only clay or reading restart file )
631      !!
632      !!   History :
633      !!        !  06-07  (C. Ethe)  original
634      !!----------------------------------------------------------------------
635 
636      ! local variables
637      INTEGER :: &
638         ji, jk, zhipor
639
640      !--------------------------------------------------------------------
641 
642
643      IF( .NOT. ln_rst_sed ) THEN
644
645          WRITE(numsed,*) ' Initilization of default values of sediment components'
646
647         ! default values for initial pore water concentrations [mol/l]
648         pwcp(:,:,:) = 0.
649         ! default value for initial solid component (fraction of dry weight dim=[0])
650         ! clay
651         solcp(:,:,:) = 0.
652         solcp(:,2:jpksed,jsclay) = 1.
653
654         ! Initialization of [h+] and [co3--]
655
656         zhipor = 8.
657         ! Initialization of [h+] in mol/kg
658         DO jk = 1, jpksed
659            DO ji = 1, jpoce
660               hipor (ji,jk) = 10.**( -1. * zhipor )
661            ENDDO
662         ENDDO
663
664         co3por(:,:) = 0.
665
666      ELSE   
667 
668         WRITE(numsed,*) ' Initilization of Sediment components from restart'
669
670         CALL sed_rst_read
671
672      ENDIF
673
674
675      ! Load initial Pisces Data for bot. wat. Chem and fluxes
676      CALL sed_dta ( nitsed000 ) 
677
678      ! Initialization of chemical constants
679      CALL sed_chem ( nitsed000 )
680
681      ! Stores initial sediment data for mass balance calculation
682      pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat ) 
683      solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol) 
684
685      ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density)
686      DO jk = 1, jpksed
687         hipor(1:jpoce,jk) = hipor(1:jpoce,jk) * densSW(1:jpoce)
688      ENDDO
689
690
691      ! In default case - no restart - sedco3 is run to initiate [h+] and [co32-]
692      ! Otherwise initiate values of pH and co3 read in restart
693      IF( .NOT. ln_rst_sed ) THEN
694         ! sedco3 is run to initiate[h+] [co32-] in mol/l of solution
695         CALL sed_co3 ( nitsed000 )
696
697      ENDIF
698           
699   END SUBROUTINE sed_init_data
700
701   SUBROUTINE sed_init_wri
702
703      INTEGER :: jk
704
705      WRITE(numsed,*)' '
706      WRITE(numsed,*)'======== Write summary of initial state   ============'
707      WRITE(numsed,*)' '
708      WRITE(numsed,*)' '
709      WRITE(numsed,*)'-------------------------------------------------------------------'
710      WRITE(numsed,*)' Initial Conditions '
711      WRITE(numsed,*)'-------------------------------------------------------------------'
712      WRITE(numsed,*)'dzm = dzkbot minimum to calculate ', 0.
713      WRITE(numsed,*)'Local zone : jpi, jpj : ',jpi, jpj
714      WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce
715      WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1)
716      WRITE(numsed,*)'Coeff diff for k=1 (cm2/s) : ',diff(1)
717      WRITE(numsed,*)' nb solid comp : ',jpsol
718      WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3)'
719      WRITE(numsed,*)'weight mol 1,2,3,4'
720      WRITE(numsed,'(4(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal)
721      WRITE(numsed,*)'nb dissolved comp',jpwat
722      WRITE(numsed,*)'(1=silicic acid,2="dissolved" clay,3=O2,4=DIC,5=Nitrate,&
723         &6=Phosphates,7=Alk))'
724      WRITE(numsed,*)'Psat (umol/l) for silicic Acid and "dissolved" clay'
725      WRITE(numsed,'(2(F0.2,3X))') sat_sil * 1e+6,  sat_clay * 1e+6
726      WRITE(numsed,*)'reaction rate rc for Op/si,Clay,POC/O2,caco3, POC/No3 (an-1)'
727      WRITE(numsed,'(5(F0.2,3X))') reac_sil * ryear, reac_clay * ryear, reac_poc * ryear, &
728                                   reac_cal * ryear, reac_no3 * ryear
729      WRITE(numsed,*)'redfield coef C,O,N P Dit '
730      WRITE(numsed,'(5(F0.2,3X))')1./spo4r,so2ut/spo4r,srno3/spo4r,spo4r/spo4r,srDnit/spo4r
731      WRITE(numsed,*)'threshold for stating denitrification [mol/l]'
732      WRITE(numsed,'(1PE8.2)') sthrO2
733      WRITE(numsed,*)'-------------------------------------------------------------------'
734      WRITE(numsed,*)'Min-Max-Mean'
735      WRITE(numsed,*)'For each variable : min, max, moy value observed on selected local zone'
736      WRITE(numsed,*)'-------------------------------------------------------------------'
737      WRITE(numsed,*)'thickness of the last wet layer dzkbot [m]'
738      WRITE(numsed,'(3(F0.2,3X))') MINVAL(dzkbot(1:jpoce))/100.,MAXVAL(dzkbot(1:jpoce))/100.,&
739         &SUM(dzkbot(1:jpoce))/jpoce/100.
740      WRITE(numsed,*)'temp [°C]'
741      WRITE(numsed,'(3(F0.2,3X))') MINVAL(temp(1:jpoce)),MAXVAL(temp(1:jpoce)),&
742         &                         SUM(temp(1:jpoce))/jpoce
743      WRITE(numsed,*)'salt o/oo'
744      WRITE(numsed,'(3(F0.2,3X))')MINVAL(salt(1:jpoce)),MAXVAL(salt(1:jpoce)),&
745         &                        SUM(salt(1:jpoce))/jpoce
746#if defined key_sed_off
747      WRITE(numsed,*)'pressure [bar] (depth in m is about 10*pressure)'
748      WRITE(numsed,'(3(F0.2,3X))') MINVAL(press(1:jpoce)),MAXVAL(press(1:jpoce)),&
749         &                         SUM(press(1:jpoce))/jpoce
750#endif
751      WRITE(numsed,*)'density of Sea Water'
752      WRITE(numsed,'(3(F0.2,3X))') MINVAL(densSW(1:jpoce)), MAXVAL(densSW(1:jpoce)),&
753         &                         SUM(densSW(1:jpoce))/jpoce
754      WRITE(numsed,*)''
755      WRITE(numsed,*)'     Dissolved Components '
756      WRITE(numsed,*)'     ====================='
757      WRITE(numsed,*)'[Si(OH)4] dissolved (1)(k=1)(µmol/l)(and min value in mol/kg  of solution)'
758      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,&
759         &                         SUM(pwcp(1:jpoce,1,jwsil))*1.e+6/jpoce,&
760         &                         MINVAL(pwcp(1:jpoce,1,jwsil)*1.e+6/densSW(1:jpoce))
761      WRITE(numsed,*)'[O2]    dissolved (3) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
762      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,&
763         &SUM(pwcp(1:jpoce,1,jwoxy))*1.e+6/jpoce,&
764         &MINVAL(pwcp(1:jpoce,1,jwoxy)*1.e+6/densSW(1:jpoce))
765      WRITE(numsed,*)'[DIC]    dissolved (4) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
766      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,&
767         &SUM(pwcp(1:jpoce,1,jwdic))*1.e+6/jpoce,&
768         &MINVAL(pwcp(1:jpoce,1,jwdic)*1.e+6/densSW(1:jpoce))
769      WRITE(numsed,*)'[NO3]    dissolved (5) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
770      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,&
771         &SUM(pwcp(1:jpoce,1,jwno3))*1.e+6/jpoce,&
772         &MINVAL(pwcp(1:jpoce,1,jwno3)*1.e+6/densSW(1:jpoce))
773      WRITE(numsed,*)'[PO4]    dissolved (6) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
774      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,&
775         &SUM(pwcp(1:jpoce,1,jwpo4))*1.e+6/jpoce,&
776         &MINVAL(pwcp(1:jpoce,1,jwpo4)*1.e+6/densSW(1:jpoce))
777      WRITE(numsed,*)'[Alk]    dissolved (7) (k=1)(µequi)(and min value in mol/kg  of solution)'
778      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,&
779         &SUM(pwcp(1:jpoce,1,jwalk))*1.e+6/jpoce,&
780         &MINVAL(pwcp(1:jpoce,1,jwalk)*1.e+6/densSW(1:jpoce))
781      WRITE(numsed,*)'[DIC13]  dissolved (8) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
782      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,&
783         &SUM(pwcp(1:jpoce,1,jwc13))*1.e+6/jpoce,&
784         &MINVAL(pwcp(1:jpoce,1,jwc13)*1.e+6/densSW(1:jpoce))
785      WRITE(numsed,*)''
786      WRITE(numsed,*)'     Solid Components '
787      WRITE(numsed,*)'     ====================='
788      WRITE(numsed,*)'nmole of Opale rained per dt'
789      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsopal))*dtsed,MAXVAL(rainrm(1:jpoce,jsopal))*dtsed,&
790         &SUM(rainrm(1:jpoce,1))*dtsed/jpoce
791      WRITE(numsed,*)'nmole of Clay rained per dt'
792      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsclay))*dtsed,MAXVAL(rainrm(1:jpoce,jsclay))*dtsed,&
793         &SUM(rainrm(1:jpoce,jsclay))*dtsed/jpoce
794      WRITE(numsed,*)'nmole of POC rained per dt'
795      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jspoc))*dtsed,MAXVAL(rainrm(1:jpoce,jspoc))*dtsed,&
796         &SUM(rainrm(1:jpoce,jspoc))*dtsed/jpoce
797      WRITE(numsed,*)'nmole of CaCO3 rained per dt'
798      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jscal))*dtsed,MAXVAL(rainrm(1:jpoce,jscal))*dtsed,&
799         &SUM(rainrm(1:jpoce,jscal))*dtsed/jpoce
800      WRITE(numsed,*)' '
801      WRITE(numsed,*)'Weight frac of opal rained (%) '
802      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
803         &MAXVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
804         & SUM(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100./jpoce
805      WRITE(numsed,*)'Weight frac of clay  rained (%) '
806      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
807         &MAXVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
808         &SUM(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100./jpoce
809      WRITE(numsed,*)'Weight frac of POC rained (%)'
810      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
811         &MAXVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
812         &SUM(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100./jpoce
813      WRITE(numsed,*)'Weight frac of CaCO3  rained (%)'
814      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
815         &MAXVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
816         &SUM(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100./jpoce
817      WRITE(numsed,*)''
818      WRITE(numsed,*)' Thickness of sediment layer added by particule rain, dzdep cm '
819      WRITE(numsed,*)' =============================================================='
820      WRITE(numsed,'(3(F0.5,2X))') MINVAL(dzdep(1:jpoce)),MAXVAL(dzdep(1:jpoce)),SUM(dzdep(1:jpoce))/jpoce
821      WRITE(numsed,*)''
822      WRITE(numsed,*)' chemical constants K1,pK1,K2,pK2,Kw,pKw and Kb pKb (min max) [mol/kgsol] or [mol/kgsol]2 '
823      WRITE(numsed,*)' ========================================================================================='
824      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak1s(1:jpoce)),MAXVAL(ak1s(1:jpoce)),-LOG10(MINVAL(ak1s(1:jpoce))),&
825         &-LOG10(MAXVAL(ak1s(1:jpoce)))
826      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak2s(1:jpoce)),MAXVAL(ak2s(1:jpoce)),-LOG10(MINVAL(ak2s(1:jpoce))),&
827         &-LOG10(MAXVAL(ak2s(1:jpoce)))
828      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akws(1:jpoce)),MAXVAL(akws(1:jpoce)),-LOG10(MINVAL(akws(1:jpoce))),&
829         &-LOG10(MAXVAL(akws(1:jpoce)))
830      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akbs(1:jpoce)),MAXVAL(akbs(1:jpoce)),-LOG10(MINVAL(akbs(1:jpoce))),&
831         &-LOG10(MAXVAL(akbs(1:jpoce)))
832      WRITE(numsed,*)'and boron'
833      WRITE(numsed,'(2(1PE10.3,2X))')MINVAL(borats(1:jpoce)),MAXVAL(borats(1:jpoce))
834      WRITE(numsed,*)''
835      WRITE(numsed,*)' Compo of initial sediment for point jpoce'
836      WRITE(numsed,*)' ========================================='
837      WRITE(numsed,*)'solcp(1),   solcp(2),   solcp(3),   solcp(4),   hipor,      pH,         co3por'
838      DO jk = 1,jpksed
839         WRITE(numsed,'(7(1PE10.3,2X))')solcp(jpoce,jk,jsopal),solcp(jpoce,jk,jsclay),solcp(jpoce,jk,jspoc),solcp(jpoce,jk,jscal),&
840            &       hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk)
841      ENDDO
842      WRITE(numsed,'(A82)')'pwcp(1),    pwcp(2),    pwcp(3),    pwcp(4),    pwcp(5),    pwcp(6),    pwcp(7)'
843      DO jk = 1, jpksed
844         WRITE(numsed,'(7(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),pwcp(jpoce,jk,jwdic),&
845            & pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13)
846      ENDDO
847      WRITE(numsed,*) ' '
848      WRITE(numsed,*) ' End Of Initialization '
849      WRITE(numsed,*) ' '
850!
851   END SUBROUTINE sed_init_wri
852#else
853   !!----------------------------------------------------------------------
854   !!   Dummy module :                      NO Sediment model
855   !!----------------------------------------------------------------------
856CONTAINS
857   SUBROUTINE sed_ini              ! Empty routine
858   END SUBROUTINE sed_ini
859#endif
860
861
862END MODULE sedini
Note: See TracBrowser for help on using the repository browser.