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

Last change on this file since 4528 was 4292, checked in by cetlod, 7 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

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