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

source: trunk/NEMOGCM/NEMO/TOP_SRC/SED/sedini.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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