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 @ 1312

Last change on this file since 1312 was 1290, checked in by cetlod, 15 years ago
open files with subroutine ctlopn in sediment model, see ticket:319
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	CHARACTER (len=40) :: csedout
79	#if defined key_sed_off
80	INTEGER :: numblt
81	INTEGER :: nummsh
82	REAL(wp) , DIMENSION(jpi,jpj) :: zdta
83	#endif
84	!!----------------------------------------------------------------------
85
86
87	! Reading namelist.sed variables
88	!---------------------------------------
89
90	csedout = 'sediment.output'
91	CALL ctlopn( numsed, csedout, 'OLD', 'FORMATTED', 'SEQUENTIAL', &
92	& 1, numout, .FALSE., 1 )
93
94	WRITE(numsed,*)
95	WRITE(numsed,*) ' L S C E - C E A'
96	WRITE(numsed,*) ' SEDIMENT model'
97	WRITE(numsed,*) ' version 2.0 (2007) '
98	WRITE(numsed,*)
99	WRITE(numsed,*)
100
101
102	WRITE(numsed,*) ' sed_init : Initialization of sediment module '
103	WRITE(numsed,*) ' '
104
105	! Determination of sediments number of points and allocate global variables
106	#if defined key_sed_off
107	! Reading Netcdf Pisces file for deepest water layer thickness [m]
108	! This data will be used as mask to merdge space in 1D vector
109	!----------------------------------------------------------------
110
111	CALL iom_open ( 'mesh_mask' , nummsh )
112	CALL iom_open ( 'e3tbot' , numblt )
113
114	! mask sediment points for outputs
115	CALL iom_get( nummsh, jpdom_data, 'tmask' , tmask )
116	CALL iom_get( nummsh, jpdom_data, 'mbathy', sbathy )
117
118	! longitude/latitude values
119	CALL iom_get ( nummsh, jpdom_data,'nav_lon', glamt(:,:) )
120	CALL iom_get ( nummsh, jpdom_data,'nav_lat', gphit(:,:) )
121
122	! bottom layer thickness
123	CALL iom_get ( numblt, jpdom_data, 'E3TBOT', zdta(:,:) )
124	epkbot(:,:) = 0.
125	DO jj = 1, jpj
126	DO ji = 1, jpi
127	ikt = MAX( INT( sbathy(ji,jj) ) - 1, 1 )
128	IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = zdta(ji,jj)
129	ENDDO
130	ENDDO
131
132	CALL iom_close( nummsh )
133	CALL iom_close( numblt )
134	#else
135
136	epkbot(:,:) = 0.
137	DO jj = 1, jpj
138	DO ji = 1, jpi
139	ikt = MAX( mbathy(ji,jj) - 1, 1 )
140	IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ikt)
141	ENDDO
142	ENDDO
143	#endif
144
145
146	! computation of total number of ocean points
147	!--------------------------------------------
148	jpoce = COUNT( epkbot(:,:) > 0. )
149
150
151	! Allocate memory size of global variables
152	ALLOCATE( pwcp (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp_dta (jpoce,jpwat) )
153	ALLOCATE( solcp(jpoce,jpksed,jpsol) ) ; ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ; ALLOCATE( rainrm_dta(jpoce,jpsol) )
154	ALLOCATE( rainrm(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) )
155	ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) )
156	ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) )
157	ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) )
158	ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) )
159	ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) )
160
161	! Initialization of global variables
162	pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0.
163	solcp (:,:,:) = 0. ; solcp0(:,:,:) = 0. ; rainrm_dta(:,:) = 0.
164	rainrm(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0.
165	dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0.
166	temp (: ) = 0. ; salt (: ) = 0.
167	press (: ) = 0. ; densSW (: ) = 0.
168	hipor (:,: ) = 0. ; co3por (:,: ) = 0.
169	dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0.
170
171	! Chemical variables
172	ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) )
173	ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) )
174	ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) )
175	ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) )
176
177	akbs (:) = 0. ; ak1s (:) = 0. ; ak2s (:) = 0. ; akws (:) = 0.
178	ak1ps (:) = 0. ; ak2ps (:) = 0. ; ak3ps (:) = 0. ; aksis (:) = 0.
179	aksps (:) = 0. ; ak12s (:) = 0. ; ak12ps(:) = 0. ; ak123ps(:) = 0.
180	borats(:) = 0. ; calcon2(:) = 0.
181
182	! Mass balance calculation
183	ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ; ALLOCATE( rloss(jpoce, jpsol) )
184	ALLOCATE( tokbot (jpoce, jpwat) )
185
186	fromsed(:,:) = 0. ; tosed(:,:) = 0. ; rloss(:,:) = 0. ; tokbot(:,:) = 0.
187
188	! Read sediment Namelist
189	!-------------------------
190	CALL sed_init_nam
191
192	! Initialization of sediment geometry
193	!------------------------------------
194	CALL sed_init_geom
195
196
197	! sets initial sediment composition
198	! ( only clay or reading restart file )
199	!---------------------------------------
200	CALL sed_init_data
201
202
203	CALL sed_init_wri
204
205
206	END SUBROUTINE sed_init
207
208
209	SUBROUTINE sed_init_geom
210	!!----------------------------------------------------------------------
211	!! * ROUTINE sed_init_geom *
212	!!
213	!! ** Purpose : Initialization of sediment geometry
214	!! - Read the deepest water layer thickness
215	!! ( using as mask ) in Netcdf file
216	!! - sets sediment grid, porosity and molecular weight
217	!! and others constants
218	!!
219	!! History :
220	!! ! 06-07 (C. Ethe) Original
221	!!----------------------------------------------------------------------
222	!! * Modules used
223	!! * local declarations
224	INTEGER :: &
225	ji, jj, jk
226
227	#if defined key_sed_off
228	REAL(wp) , DIMENSION (jpi,jpj) :: zdta
229	INTEGER :: numpres
230	#endif
231	!----------------------------------------------------------
232
233	WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
234	WRITE(numsed,*) ' '
235
236	! Computation of 1D array of sediments points
237	indoce = 0
238	DO jj = 1, jpj
239	DO ji = 1, jpi
240	IF ( epkbot(ji,jj) > 0. ) THEN
241	indoce = indoce + 1
242	iarroce(indoce) = (jj - 1) * jpi + ji
243	ENDIF
244	END DO
245	END DO
246
247	IF( indoce .NE. jpoce ) THEN
248	WRITE(numsed,*) ' '
249	WRITE(numsed,*) 'Warning : number of ocean points indoce = ', indoce, &
250	& ' doesn''t match number of point where epkbot>0 jpoce = ', jpoce
251	WRITE(numsed,*) ' '
252	WRITE(numsed,*) ' '
253	STOP
254	ELSE
255	WRITE(numsed,*) ' '
256	WRITE(numsed,*) ' total number of ocean points jpoce = ',jpoce
257	WRITE(numsed,*) ' '
258	ENDIF
259
260	#if defined key_sed_off
261
262	! Reading Netcdf Pisces file for deepest water layer thickness [m]
263	! This data will be used as mask to merdge space in 1D vector
264	!----------------------------------------------------------------
265	CALL iom_open ( 'pressbot' , numpres )
266
267	! pressure in bars
268	CALL iom_get ( numpres, jpdom_data,'BATH', zdta(:,:) )
269	CALL pack_arr ( jpoce, press(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) )
270	press(1:jpoce) = press(1:jpoce) * 1.025e-1
271
272	CALL iom_close ( numpres )
273	#endif
274
275
276	! mask sediment points for outputs
277	DO jk = 1, jpksed
278	tmasksed(:,:,jk) = tmask(:,:,1)
279	ENDDO
280
281	! initialization of dzkbot in [cm]
282	!------------------------------------------------
283	CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
284	dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2
285
286	! Geometry and constants
287	! sediment layer thickness [cm]
288	! (1st layer= diffusive layer = pur water)
289	!------------------------------------------
290	dz(1) = 0.1
291	dz(2) = 0.3
292	dz(3) = 0.3
293	dz(4) = 0.5
294	dz(5) = 0.5
295	dz(6) = 0.5
296	dz(7) = 1.
297	dz(8) = 1.
298	dz(9) = 1.
299	dz(10) = 2.45
300	dz(11) = 2.45
301
302	DO jk = 1, jpksed
303	DO ji = 1, jpoce
304	dz3d(ji,jk) = dz(jk)
305	END DO
306	ENDDO
307
308	! Depth(jk)= depth of middle of each layer
309	!----------------------------------------
310	profsed(1) = -dz(1)/ 2.
311	DO jk = 2, jpksed
312	profsed(jk) = profsed(jk-1) + dz(jk-1) / 2. + dz(jk) / 2.
313	ENDDO
314
315	! Porosity profile [0]
316	!---------------------
317	por(1) = 1.
318	por(2) = 0.95
319	por(3) = 0.9
320	por(4) = 0.85
321	por(5) = 0.81
322	por(6) = 0.75
323	por(7) = 0.75
324	por(8) = 0.75
325	por(9) = 0.75
326	por(10) = 0.75
327	por(11) = 0.75
328
329	! inverse of Porosity profile
330	!-----------------------------
331	por1(:) = 1. - por(:)
332
333	! Volumes of pore water and solid fractions (vector and array)
334	! WARNING : volw(1) and vols(1) are sublayer volums
335	volw(:) = dz(:) * por(:)
336	vols(:) = dz(:) * por1(:)
337
338	! temporary new value for dz3d(:,1)
339	dz3d(1:jpoce,1) = dzkbot(1:jpoce)
340
341	! WARNING : volw3d(:,1) and vols3d(:,1) are deepest water column volums
342	DO jk = 1, jpksed
343	volw3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por (jk)
344	vols3d(1:jpoce,jk) = dz3d(1:jpoce,jk) * por1(jk)
345	ENDDO
346
347	! Back to the old sublayer vlaue for dz3d(:,1)
348	dz3d(1:jpoce,1) = dz(1)
349
350
351	!---------------------------------------------
352	! Molecular weight [g/mol] for solid species
353	!---------------------------------------------
354
355
356	! opal=sio20.4(h20)=28+216+0.4*(2+16)
357	!---------------------------------------
358	mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. )
359
360	! clay
361	! some kind of Illit (according to Pape)
362	! K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2\|(Si3.41Al0.59)O10]
363	!--------------------------------------------------------------------
364	mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &
365	& 0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. + &
366	& 0.59 * 27. + 10. * 16.
367
368	! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
369	! But den sity of Poc is an Hydrated material (= POC + 30H2O)
370	! So C(122)H(355)O(120)N(16)P(1)
371	!------------------------------------------------------------
372	mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
373	& 16. * 14. + 31. ) / 122.
374
375	! CaCO3
376	!---------
377	mol_wgt(jscal) = 40. + 12. + 3. * 16.
378
379	! Density of solid material in sediment [g/cm**3]
380	!------------------------------------------------
381	dens = 2.6
382
383	! Initialization of diffusion coefficient as function of porosity [cm**2/s]
384	!--------------------------------------------------------------------
385	diff(:) = dcoef * por(:)
386
387
388	! Initialization of time step as function of porosity [cm**2/s]
389	!------------------------------------------------------------------
390	rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) )
391
392	END SUBROUTINE sed_init_geom
393
394	SUBROUTINE sed_init_nam
395	!!----------------------------------------------------------------------
396	!! * ROUTINE sed_init_nam *
397	!!
398	!! ** Purpose : Initialization of sediment geometry
399	!! - Reading namelist and defines constants variables
400	!!
401	!! History :
402	!! ! 06-07 (C. Ethe) Original
403	!!----------------------------------------------------------------------
404
405	INTEGER :: &
406	numnamsed = 28
407
408	TYPE PSED
409	CHARACTER(len = 20) :: snamesed !: short name
410	CHARACTER(len = 80 ) :: lnamesed !: long name
411	CHARACTER(len = 20 ) :: unitsed !: unit
412	END TYPE PSED
413
414	TYPE(PSED) , DIMENSION(jpsol ) :: sedsol
415	TYPE(PSED) , DIMENSION(jpwat ) :: sedwat
416	TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
417	TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
418
419	NAMELIST/nam_time/nfreq
420	NAMELIST/nam_trased/sedsol, sedwat
421	NAMELIST/nam_diased/seddiag3d, seddiag2d
422	NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef
423	NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, &
424	& rcorg, o2seuil, rcorgN
425	NAMELIST/nam_cal/rccal
426	NAMELIST/nam_dc13/pdb, rc13P, rc13Ca
427	NAMELIST/nam_btb/dbiot
428	NAMELIST/nam_rst/ln_rst_sed
429
430	INTEGER :: jn, jn1
431	CHARACTER(len=40) :: csednam
432	!-------------------------------------------------------
433
434	WRITE(numsed,*) ' sed_init_nam : Read namelists '
435	WRITE(numsed,*) ' '
436
437	! ryear = 1 year converted in second
438	!------------------------------------
439	WRITE(numsed,*) ' '
440	WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
441	WRITE(numsed,*) ' '
442
443	! Reading namelist.sed variables
444	!---------------------------------
445	csednam = 'namelist.sediment'
446	CALL ctlopn( numnamsed, csednam, 'OLD', 'FORMATTED', 'SEQUENTIAL', &
447	& 1, numout, .FALSE., 1 )
448
449	dtsed = rdt
450	#if ! defined key_sed_off
451	nitsed000 = nittrc000
452	nitsedend = nitend
453	nwrised = nwritetrc
454	#else
455	nitsed000 = nit000
456	nitsedend = nitend
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 = rc13PPDB rc13P = ', rc13P
602	WRITE(numsed,) ' 13C/12C in CaCO3 = rc13CaPDB 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 10pressure)'
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	CONTAINS
860	SUBROUTINE sed_ini ! Empty routine
861	END SUBROUTINE sed_ini
862	#endif
863
864
865	END MODULE sedini

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

Download in other formats: