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/NEMOGCM/NEMO/TOP_SRC/SED – NEMO

Context Navigation

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

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 34.1 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) , DIMENSION(jpi,jpj) :: 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	! 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 )
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 = mbkt(ji,jj)
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
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	TYPE PSED
406	CHARACTER(len = 20) :: snamesed !: short name
407	CHARACTER(len = 80 ) :: lnamesed !: long name
408	CHARACTER(len = 20 ) :: unitsed !: unit
409	END TYPE PSED
410
411	TYPE(PSED) , DIMENSION(jpsol ) :: sedsol
412	TYPE(PSED) , DIMENSION(jpwat ) :: sedwat
413	TYPE(PSED) , DIMENSION(jpdia3dsed) :: seddiag3d
414	TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
415
416	NAMELIST/nam_time/nfreq
417	NAMELIST/nam_trased/sedsol, sedwat
418	NAMELIST/nam_diased/seddiag3d, seddiag2d
419	NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef
420	NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, &
421	& rcorg, o2seuil, rcorgN
422	NAMELIST/nam_cal/rccal
423	NAMELIST/nam_dc13/pdb, rc13P, rc13Ca
424	NAMELIST/nam_btb/dbiot
425	NAMELIST/nam_rst/ln_rst_sed
426
427	INTEGER :: jn, jn1
428	!-------------------------------------------------------
429
430	WRITE(numsed,*) ' sed_init_nam : Read namelists '
431	WRITE(numsed,*) ' '
432
433	! ryear = 1 year converted in second
434	!------------------------------------
435	WRITE(numsed,*) ' '
436	WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
437	WRITE(numsed,*) ' '
438
439	! Reading namelist.sed variables
440	!---------------------------------
441	CALL ctl_opn( numnamsed, 'namelist.sediment', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
442
443	dtsed = rdt
444	nitsed000 = nit000
445	nitsedend = nitend
446	#if ! defined key_sed_off
447	nwrised = nwritetrc
448	#else
449	nwrised = nwrite
450	#endif
451	! Diffraction/reaction parameters
452	!----------------------------------
453	READ( numnamsed, nam_time )
454	WRITE(numsed,*) ' namelist nam_time'
455
456	#if ! defined key_sed_off
457	nfreq = 1
458	#endif
459
460	WRITE(numsed,*) ' sedimentation time step dtsed = ', dtsed
461	WRITE(numsed,*) ' 1st time step for sediment. nitsed000 = ', nitsed000
462	WRITE(numsed,*) ' last time step for sediment. nitsedend = ', nitsedend
463	WRITE(numsed,*) ' frequency of sediment outputs nwrised = ', nwrised
464	WRITE(numsed,*) ' frequency of restoring inputs data nfreq = ', nfreq
465	WRITE(numsed,*) ' '
466
467	REWIND( numnamsed ) ! read nattrc
468	READ ( numnamsed, nam_trased )
469
470	DO jn = 1, jpsol
471	sedtrcd(jn) = sedsol(jn)%snamesed
472	sedtrcl(jn) = sedsol(jn)%lnamesed
473	sedtrcu(jn) = sedsol(jn)%unitsed
474	END DO
475
476	DO jn = 1, jpwat
477	jn1 = jn + jpsol
478	sedtrcd(jn1) = sedwat(jn)%snamesed
479	sedtrcl(jn1) = sedwat(jn)%lnamesed
480	sedtrcu(jn1) = sedwat(jn)%unitsed
481	END DO
482
483	WRITE(numsed,*) ' namelist nam_trased'
484	WRITE(numsed,*) ' '
485	DO jn = 1, jptrased
486	WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
487	WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
488	WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
489	WRITE(numsed,*) ' '
490	END DO
491	WRITE(numsed,*) ' '
492
493
494	REWIND( numnamsed )
495	READ( numnamsed, nam_diased )
496
497	DO jn = 1, jpdia3dsed
498	seddia3d(jn) = seddiag3d(jn)%snamesed
499	seddia3l(jn) = seddiag3d(jn)%lnamesed
500	seddia3u(jn) = seddiag3d(jn)%unitsed
501	END DO
502
503	DO jn = 1, jpdia2dsed
504	seddia2d(jn) = seddiag2d(jn)%snamesed
505	seddia2l(jn) = seddiag2d(jn)%lnamesed
506	seddia2u(jn) = seddiag2d(jn)%unitsed
507	END DO
508
509	WRITE(numsed,*) ' namelist nam_diased'
510	WRITE(numsed,*) ' '
511	DO jn = 1, jpdia3dsed
512	WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
513	WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
514	WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
515	WRITE(numsed,*) ' '
516	END DO
517
518	DO jn = 1, jpdia2dsed
519	WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
520	WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
521	WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
522	WRITE(numsed,*) ' '
523	END DO
524
525	WRITE(numsed,*) ' '
526
527
528	! Diffraction/reaction parameters
529	!----------------------------------
530	REWIND( numnamsed )
531	READ( numnamsed, nam_reac )
532	WRITE(numsed,*) ' namelist nam_reac'
533	WRITE(numsed,*) ' saturation for si sisat = ', sisat
534	WRITE(numsed,*) ' saturation for clay claysat = ', claysat
535	WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal
536	WRITE(numsed,*) ' reactivity for clay rcclay = ', rcclay
537	WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef
538	WRITE(numsed,*) ' '
539
540
541	! Unity conversion to get saturation conc. psat in [mol.l-1]
542	! and reactivity rc in [l.mol-1.s-1]
543	!----------------------------------------------------------
544	sat_sil = sisat * 1.e-6
545	sat_clay = claysat * 1.e-6
546
547	reac_sil = rcopal / ryear
548	reac_clay = rcclay / ryear
549
550
551	! Additional parameter linked to POC/O2/No3/Po4
552	!----------------------------------------------
553	REWIND( numnamsed )
554	READ( numnamsed, nam_poc )
555	WRITE(numsed,*) ' namelist nam_poc'
556	WRITE(numsed,*) ' Redfield coef for oxy redO2 = ', redO2
557	WRITE(numsed,*) ' Redfield coef for no3 redNo3 = ', redNo3
558	WRITE(numsed,*) ' Redfield coef for po4 redPo4 = ', redPo4
559	WRITE(numsed,*) ' Redfield coef for carbone redC = ', redC
560	WRITE(numsed,*) ' Redfield coef for denitri redDnit = ', redDnit
561	WRITE(numsed,*) ' reactivity for POC/O2 rcorg = ', rcorg
562	WRITE(numsed,*) ' threshold O2 concen. o2seuil = ', o2seuil
563	WRITE(numsed,*) ' reactivity for POC/NO3 rcorgN = ', rcorgN
564	WRITE(numsed,*) ' '
565
566
567	so2ut = redO2 / redC
568	srno3 = redNo3 / redC
569	spo4r = redPo4 / redC
570	srDnit = redDnit / redC
571	sthro2 = o2seuil * 1.e-6 ! threshold O2 concen. in [mol.l-1]
572	! reactivity rc in [l.mol-1.s-1]
573	reac_poc = rcorg / ryear
574	reac_no3 = rcorgN / ryear
575
576
577	! Carbonate parameters
578	!---------------------
579	READ( numnamsed, nam_cal )
580	WRITE(numsed,*) ' namelist nam_cal'
581	WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal
582	WRITE(numsed,*) ' '
583
584	! reactivity rc in [l.mol-1.s-1]
585	reac_cal = rccal / ryear
586
587
588	! C13 parameters
589	!----------------
590	READ( numnamsed, nam_dc13 )
591	WRITE(numsed,*) ' namelist nam_dc13 '
592	WRITE(numsed,*) ' 13C/12C in PD Belemnite PDB = ', pdb
593	WRITE(numsed,) ' 13C/12C in POC = rc13PPDB rc13P = ', rc13P
594	WRITE(numsed,) ' 13C/12C in CaCO3 = rc13CaPDB rc13Ca = ', rc13Ca
595	WRITE(numsed,*) ' '
596
597
598	! Bioturbation parameter
599	!------------------------
600	READ( numnamsed, nam_btb )
601	WRITE(numsed,*) ' namelist nam_btb '
602	WRITE(numsed,*) ' coefficient for bioturbation dbiot = ', dbiot
603	WRITE(numsed,*) ' '
604
605	! Unity convertion to get bioturb coefficient in [cm2.s-1]
606	db = dbiot / ( ryear * 1000. )
607
608	! Initial value (t=0) for sediment pore water and solid components
609	!----------------------------------------------------------------
610	READ( numnamsed, nam_rst )
611	WRITE(numsed,*) ' namelist nam_rst '
612	WRITE(numsed,*) ' boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed
613	WRITE(numsed,*) ' '
614
615	CLOSE( numnamsed )
616
617	END SUBROUTINE sed_init_nam
618
619	SUBROUTINE sed_init_data
620	!!----------------------------------------------------------------------
621	!! * ROUTINE sed_init_data *
622	!!
623	!! ** Purpose : Initialization of sediment module
624	!! - sets initial sediment composition
625	!! ( only clay or reading restart file )
626	!!
627	!! History :
628	!! ! 06-07 (C. Ethe) original
629	!!----------------------------------------------------------------------
630
631	! local variables
632	INTEGER :: &
633	ji, jk, zhipor
634
635	!--------------------------------------------------------------------
636
637
638	IF( .NOT. ln_rst_sed ) THEN
639
640	WRITE(numsed,*) ' Initilization of default values of sediment components'
641
642	! default values for initial pore water concentrations [mol/l]
643	pwcp(:,:,:) = 0.
644	! default value for initial solid component (fraction of dry weight dim=[0])
645	! clay
646	solcp(:,:,:) = 0.
647	solcp(:,2:jpksed,jsclay) = 1.
648
649	! Initialization of [h+] and [co3--]
650
651	zhipor = 8.
652	! Initialization of [h+] in mol/kg
653	DO jk = 1, jpksed
654	DO ji = 1, jpoce
655	hipor (ji,jk) = 10.*( -1. zhipor )
656	ENDDO
657	ENDDO
658
659	co3por(:,:) = 0.
660
661	ELSE
662
663	WRITE(numsed,*) ' Initilization of Sediment components from restart'
664
665	CALL sed_rst_read
666
667	ENDIF
668
669
670	! Load initial Pisces Data for bot. wat. Chem and fluxes
671	CALL sed_dta ( nitsed000 )
672
673	! Initialization of chemical constants
674	CALL sed_chem ( nitsed000 )
675
676	! Stores initial sediment data for mass balance calculation
677	pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat )
678	solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol)
679
680	! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density)
681	DO jk = 1, jpksed
682	hipor(1:jpoce,jk) = hipor(1:jpoce,jk) * densSW(1:jpoce)
683	ENDDO
684
685
686	! In default case - no restart - sedco3 is run to initiate [h+] and [co32-]
687	! Otherwise initiate values of pH and co3 read in restart
688	IF( .NOT. ln_rst_sed ) THEN
689	! sedco3 is run to initiate[h+] [co32-] in mol/l of solution
690	CALL sed_co3 ( nitsed000 )
691
692	ENDIF
693
694	END SUBROUTINE sed_init_data
695
696	SUBROUTINE sed_init_wri
697
698	INTEGER :: jk
699
700	WRITE(numsed,*)' '
701	WRITE(numsed,*)'======== Write summary of initial state ============'
702	WRITE(numsed,*)' '
703	WRITE(numsed,*)' '
704	WRITE(numsed,*)'-------------------------------------------------------------------'
705	WRITE(numsed,*)' Initial Conditions '
706	WRITE(numsed,*)'-------------------------------------------------------------------'
707	WRITE(numsed,*)'dzm = dzkbot minimum to calculate ', 0.
708	WRITE(numsed,*)'Local zone : jpi, jpj : ',jpi, jpj
709	WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce
710	WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1)
711	WRITE(numsed,*)'Coeff diff for k=1 (cm2/s) : ',diff(1)
712	WRITE(numsed,*)' nb solid comp : ',jpsol
713	WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3)'
714	WRITE(numsed,*)'weight mol 1,2,3,4'
715	WRITE(numsed,'(4(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal)
716	WRITE(numsed,*)'nb dissolved comp',jpwat
717	WRITE(numsed,*)'(1=silicic acid,2="dissolved" clay,3=O2,4=DIC,5=Nitrate,&
718	&6=Phosphates,7=Alk))'
719	WRITE(numsed,*)'Psat (umol/l) for silicic Acid and "dissolved" clay'
720	WRITE(numsed,'(2(F0.2,3X))') sat_sil * 1e+6, sat_clay * 1e+6
721	WRITE(numsed,*)'reaction rate rc for Op/si,Clay,POC/O2,caco3, POC/No3 (an-1)'
722	WRITE(numsed,'(5(F0.2,3X))') reac_sil * ryear, reac_clay * ryear, reac_poc * ryear, &
723	reac_cal * ryear, reac_no3 * ryear
724	WRITE(numsed,*)'redfield coef C,O,N P Dit '
725	WRITE(numsed,'(5(F0.2,3X))')1./spo4r,so2ut/spo4r,srno3/spo4r,spo4r/spo4r,srDnit/spo4r
726	WRITE(numsed,*)'threshold for stating denitrification [mol/l]'
727	WRITE(numsed,'(1PE8.2)') sthrO2
728	WRITE(numsed,*)'-------------------------------------------------------------------'
729	WRITE(numsed,*)'Min-Max-Mean'
730	WRITE(numsed,*)'For each variable : min, max, moy value observed on selected local zone'
731	WRITE(numsed,*)'-------------------------------------------------------------------'
732	WRITE(numsed,*)'thickness of the last wet layer dzkbot [m]'
733	WRITE(numsed,'(3(F0.2,3X))') MINVAL(dzkbot(1:jpoce))/100.,MAXVAL(dzkbot(1:jpoce))/100.,&
734	&SUM(dzkbot(1:jpoce))/jpoce/100.
735	WRITE(numsed,*)'temp [°C]'
736	WRITE(numsed,'(3(F0.2,3X))') MINVAL(temp(1:jpoce)),MAXVAL(temp(1:jpoce)),&
737	& SUM(temp(1:jpoce))/jpoce
738	WRITE(numsed,*)'salt o/oo'
739	WRITE(numsed,'(3(F0.2,3X))')MINVAL(salt(1:jpoce)),MAXVAL(salt(1:jpoce)),&
740	& SUM(salt(1:jpoce))/jpoce
741	#if defined key_sed_off
742	WRITE(numsed,)'pressure [bar] (depth in m is about 10pressure)'
743	WRITE(numsed,'(3(F0.2,3X))') MINVAL(press(1:jpoce)),MAXVAL(press(1:jpoce)),&
744	& SUM(press(1:jpoce))/jpoce
745	#endif
746	WRITE(numsed,*)'density of Sea Water'
747	WRITE(numsed,'(3(F0.2,3X))') MINVAL(densSW(1:jpoce)), MAXVAL(densSW(1:jpoce)),&
748	& SUM(densSW(1:jpoce))/jpoce
749	WRITE(numsed,*)''
750	WRITE(numsed,*)' Dissolved Components '
751	WRITE(numsed,*)' ====================='
752	WRITE(numsed,*)'[Si(OH)4] dissolved (1)(k=1)(µmol/l)(and min value in mol/kg of solution)'
753	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwsil))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwsil))1.e+6,&
754	& SUM(pwcp(1:jpoce,1,jwsil))*1.e+6/jpoce,&
755	& MINVAL(pwcp(1:jpoce,1,jwsil)*1.e+6/densSW(1:jpoce))
756	WRITE(numsed,*)'[O2] dissolved (3) (k=1)(µmol/l)(and min value in mol/kg of solution)'
757	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwoxy))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwoxy))1.e+6,&
758	&SUM(pwcp(1:jpoce,1,jwoxy))*1.e+6/jpoce,&
759	&MINVAL(pwcp(1:jpoce,1,jwoxy)*1.e+6/densSW(1:jpoce))
760	WRITE(numsed,*)'[DIC] dissolved (4) (k=1)(µmol/l)(and min value in mol/kg of solution)'
761	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwdic))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwdic))1.e+6,&
762	&SUM(pwcp(1:jpoce,1,jwdic))*1.e+6/jpoce,&
763	&MINVAL(pwcp(1:jpoce,1,jwdic)*1.e+6/densSW(1:jpoce))
764	WRITE(numsed,*)'[NO3] dissolved (5) (k=1)(µmol/l)(and min value in mol/kg of solution)'
765	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwno3))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwno3))1.e+6,&
766	&SUM(pwcp(1:jpoce,1,jwno3))*1.e+6/jpoce,&
767	&MINVAL(pwcp(1:jpoce,1,jwno3)*1.e+6/densSW(1:jpoce))
768	WRITE(numsed,*)'[PO4] dissolved (6) (k=1)(µmol/l)(and min value in mol/kg of solution)'
769	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwpo4))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwpo4))1.e+6,&
770	&SUM(pwcp(1:jpoce,1,jwpo4))*1.e+6/jpoce,&
771	&MINVAL(pwcp(1:jpoce,1,jwpo4)*1.e+6/densSW(1:jpoce))
772	WRITE(numsed,*)'[Alk] dissolved (7) (k=1)(µequi)(and min value in mol/kg of solution)'
773	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwalk))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwalk))1.e+6,&
774	&SUM(pwcp(1:jpoce,1,jwalk))*1.e+6/jpoce,&
775	&MINVAL(pwcp(1:jpoce,1,jwalk)*1.e+6/densSW(1:jpoce))
776	WRITE(numsed,*)'[DIC13] dissolved (8) (k=1)(µmol/l)(and min value in mol/kg of solution)'
777	WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwc13))1.e+6,MAXVAL(pwcp(1:jpoce,1,jwc13))1.e+6,&
778	&SUM(pwcp(1:jpoce,1,jwc13))*1.e+6/jpoce,&
779	&MINVAL(pwcp(1:jpoce,1,jwc13)*1.e+6/densSW(1:jpoce))
780	WRITE(numsed,*)''
781	WRITE(numsed,*)' Solid Components '
782	WRITE(numsed,*)' ====================='
783	WRITE(numsed,*)'nmole of Opale rained per dt'
784	WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsopal))dtsed,MAXVAL(rainrm(1:jpoce,jsopal))dtsed,&
785	&SUM(rainrm(1:jpoce,1))*dtsed/jpoce
786	WRITE(numsed,*)'nmole of Clay rained per dt'
787	WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsclay))dtsed,MAXVAL(rainrm(1:jpoce,jsclay))dtsed,&
788	&SUM(rainrm(1:jpoce,jsclay))*dtsed/jpoce
789	WRITE(numsed,*)'nmole of POC rained per dt'
790	WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jspoc))dtsed,MAXVAL(rainrm(1:jpoce,jspoc))dtsed,&
791	&SUM(rainrm(1:jpoce,jspoc))*dtsed/jpoce
792	WRITE(numsed,*)'nmole of CaCO3 rained per dt'
793	WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jscal))dtsed,MAXVAL(rainrm(1:jpoce,jscal))dtsed,&
794	&SUM(rainrm(1:jpoce,jscal))*dtsed/jpoce
795	WRITE(numsed,*)' '
796	WRITE(numsed,*)'Weight frac of opal rained (%) '
797	WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
798	&MAXVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
799	& SUM(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100./jpoce
800	WRITE(numsed,*)'Weight frac of clay rained (%) '
801	WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
802	&MAXVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
803	&SUM(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100./jpoce
804	WRITE(numsed,*)'Weight frac of POC rained (%)'
805	WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
806	&MAXVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
807	&SUM(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100./jpoce
808	WRITE(numsed,*)'Weight frac of CaCO3 rained (%)'
809	WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
810	&MAXVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
811	&SUM(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100./jpoce
812	WRITE(numsed,*)''
813	WRITE(numsed,*)' Thickness of sediment layer added by particule rain, dzdep cm '
814	WRITE(numsed,*)' =============================================================='
815	WRITE(numsed,'(3(F0.5,2X))') MINVAL(dzdep(1:jpoce)),MAXVAL(dzdep(1:jpoce)),SUM(dzdep(1:jpoce))/jpoce
816	WRITE(numsed,*)''
817	WRITE(numsed,*)' chemical constants K1,pK1,K2,pK2,Kw,pKw and Kb pKb (min max) [mol/kgsol] or [mol/kgsol]2 '
818	WRITE(numsed,*)' ========================================================================================='
819	WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak1s(1:jpoce)),MAXVAL(ak1s(1:jpoce)),-LOG10(MINVAL(ak1s(1:jpoce))),&
820	&-LOG10(MAXVAL(ak1s(1:jpoce)))
821	WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak2s(1:jpoce)),MAXVAL(ak2s(1:jpoce)),-LOG10(MINVAL(ak2s(1:jpoce))),&
822	&-LOG10(MAXVAL(ak2s(1:jpoce)))
823	WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akws(1:jpoce)),MAXVAL(akws(1:jpoce)),-LOG10(MINVAL(akws(1:jpoce))),&
824	&-LOG10(MAXVAL(akws(1:jpoce)))
825	WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akbs(1:jpoce)),MAXVAL(akbs(1:jpoce)),-LOG10(MINVAL(akbs(1:jpoce))),&
826	&-LOG10(MAXVAL(akbs(1:jpoce)))
827	WRITE(numsed,*)'and boron'
828	WRITE(numsed,'(2(1PE10.3,2X))')MINVAL(borats(1:jpoce)),MAXVAL(borats(1:jpoce))
829	WRITE(numsed,*)''
830	WRITE(numsed,*)' Compo of initial sediment for point jpoce'
831	WRITE(numsed,*)' ========================================='
832	WRITE(numsed,*)'solcp(1), solcp(2), solcp(3), solcp(4), hipor, pH, co3por'
833	DO jk = 1,jpksed
834	WRITE(numsed,'(7(1PE10.3,2X))')solcp(jpoce,jk,jsopal),solcp(jpoce,jk,jsclay),solcp(jpoce,jk,jspoc),solcp(jpoce,jk,jscal),&
835	& hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk)
836	ENDDO
837	WRITE(numsed,'(A82)')'pwcp(1), pwcp(2), pwcp(3), pwcp(4), pwcp(5), pwcp(6), pwcp(7)'
838	DO jk = 1, jpksed
839	WRITE(numsed,'(7(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),pwcp(jpoce,jk,jwdic),&
840	& pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13)
841	ENDDO
842	WRITE(numsed,*) ' '
843	WRITE(numsed,*) ' End Of Initialization '
844	WRITE(numsed,*) ' '
845	!
846	END SUBROUTINE sed_init_wri
847	#else
848	!!----------------------------------------------------------------------
849	!! Dummy module : NO Sediment model
850	!!----------------------------------------------------------------------
851	CONTAINS
852	SUBROUTINE sed_ini ! Empty routine
853	END SUBROUTINE sed_ini
854	#endif
855
856
857	END MODULE sedini

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

Download in other formats: