New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sedini.F90 in trunk/NEMO/TOP_SRC/SED – NEMO

Context Navigation

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

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

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

Download in other formats: