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 branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/SED – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedini.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 10 years ago
dev_MERGE_2013 : 1st step of the merge, see ticket #1185
File size: 34.3 KB

Line
1	MODULE 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	CONTAINS
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=sio20.4(h20)=28+216+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 = rc13PPDB rc13P = ', rc13P
601	WRITE(numsed,) ' 13C/12C in CaCO3 = rc13CaPDB 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 10pressure)'
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	!!----------------------------------------------------------------------
858	CONTAINS
859	SUBROUTINE sed_ini ! Empty routine
860	END SUBROUTINE sed_ini
861	#endif
862
863
864	END MODULE sedini

Note: See TracBrowser for help on using the repository browser.

Download in other formats: