source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedini.F90 @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 5 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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