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.

p5zpoc.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z – NEMO

Context Navigation

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90 @ 7617

Last change on this file since 7617 was 7617, checked in by aumont, 8 years ago
update diagnostics + changes in quota code
Property svn:executable set to ``*
File size: 32.5 KB

Line
1	MODULE p5zpoc
2	!!======================================================================
3	!! * MODULE p5zpoc *
4	!! TOP : PISCES Compute remineralization of organic particles
5	!!=========================================================================
6	!! History : 1.0 ! 2004 (O. Aumont) Original code
7	!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
8	!! 3.4 ! 2011-06 (O. Aumont, C. Ethe) Quota model for iron
9	!! 3.6 ! 2015-05 (O. Aumont) PISCES quota
10	!! 3.6 ! 2016-03 (O. Aumont) PISCES Quota and reorganization
11	!!----------------------------------------------------------------------
12	#if defined key_pisces_quota
13	!!----------------------------------------------------------------------
14	!! 'key_top' and TOP models
15	!! 'key_pisces_quota' PISCES bio-model with variable stoichiometry
16	!!----------------------------------------------------------------------
17	!! p5z_poc : Compute remineralization of organic particles
18	!! p5z_poc_init : Initialisation of parameters for remineralisation
19	!! p5z_poc_alloc : Allocate remineralisation variables
20	!!----------------------------------------------------------------------
21	USE oce_trc ! shared variables between ocean and passive tracers
22	USE trc ! passive tracers common variables
23	USE sms_pisces ! PISCES Source Minus Sink variables
24	USE p5zsink ! Sedimentation of organic particles
25	USE prtctl_trc ! print control for debugging
26	USE iom ! I/O manager
27
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC p5z_poc ! called in p4zbio.F90
33	PUBLIC p5z_poc_init ! called in trcsms_pisces.F90
34	PUBLIC alngam
35	PUBLIC gamain
36
37	!! * Shared module variables
38	REAL(wp), PUBLIC :: xremipc !: remineralisation rate of DOC
39	REAL(wp), PUBLIC :: xremipn !: remineralisation rate of DON
40	REAL(wp), PUBLIC :: xremipp !: remineralisation rate of DOP
41	INTEGER , PUBLIC :: jcpoc !: number of lability classes
42	REAL(wp), PUBLIC :: rshape !: shape factor of the gamma distribution
43
44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: alphan, reminp
45	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: alphap
46
47	!!* Substitution
48	# include "top_substitute.h90"
49	!!----------------------------------------------------------------------
50	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
51	!! $Id: p4zrem.F90 3160 2011-11-20 14:27:18Z cetlod $
52	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
53	!!----------------------------------------------------------------------
54	CONTAINS
55
56	SUBROUTINE p5z_poc( kt, knt )
57	!!---------------------------------------------------------------------
58	!! * ROUTINE p5z_poc *
59	!!
60	!! ** Purpose : Compute remineralization of organic particles
61	!!
62	!! ** Method : - ???
63	!!---------------------------------------------------------------------
64	!
65	INTEGER, INTENT(in) :: kt, knt ! ocean time step
66	!
67	INTEGER :: ji, jj, jk, jn
68	REAL(wp) :: zremip, zremig, zdep, zstep
69	REAL(wp) :: zopon, zopop, zopoc, zopoc2, zopon2, zopop2, zofer
70	REAL(wp) :: zsizek, zsizek1, alphat, remint
71	REAL(wp) :: solgoc, zpoc
72	#if ! defined key_kriest
73	REAL(wp) :: zofer2, zofer3
74	#endif
75	REAL(wp) :: step, zrfact2
76	CHARACTER (len=25) :: charout
77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwork1
78	REAL(wp), POINTER, DIMENSION(:,: ) :: totprod, totthick, totcons
79	REAL(wp), POINTER, DIMENSION(:,:,:) :: zremipoc, zremigoc, zorem3
80	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: alphag
81	REAL(wp), POINTER, DIMENSION(:,:,:) :: zfolimi
82
83	!!---------------------------------------------------------------------
84	!
85	IF( nn_timing == 1 ) CALL timing_start('p5z_poc')
86	!
87	! Allocate temporary workspace
88	CALL wrk_alloc( jpi, jpj, jpk, zwork1 )
89	CALL wrk_alloc( jpi, jpj, totprod, totthick, totcons )
90	CALL wrk_alloc( jpi, jpj, jpk, zremipoc, zremigoc, zorem3 )
91	CALL wrk_alloc( jpi, jpj, jpk, zfolimi )
92	ALLOCATE( alphag(jpi,jpj,jpk,jcpoc) )
93
94	! Here we compute the GOC -> POC rate due to the shrinking
95	! of the fecal pellets/aggregates as a result of bacterial
96	! solubilization
97	! This is based on a fractal dimension of 2.56 and a spectral
98	! slope of -3.6 (identical to what is used in p4zsink to compute
99	! aggregation
100	solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
101
102	! Initialisation of temprary arrys
103	zremipoc(:,:,:) = xremipc
104	zremigoc(:,:,:) = xremipc
105	zorem3(:,:,:) = 0.
106	orem (:,:,:) = 0.
107
108	DO jn = 1, jcpoc
109	alphag(:,:,:,jn) = alphan(jn)
110	alphap(:,:,:,jn) = alphan(jn)
111	END DO
112
113	#if ! defined key_kriest
114	! -----------------------------------------------------------------------
115	! Lability parameterization. This is the big particles part (GOC)
116	! This lability parameterization can be activated only with the standard
117	! particle scheme. Does not work with Kriest parameterization.
118	! -----------------------------------------------------------------------
119	DO jk = 2, jpkm1
120	DO jj = 1, jpj
121	DO ji = 1, jpi
122	IF (tmask(ji,jj,jk) == 1.) THEN
123	zdep = hmld(ji,jj)
124	!
125	! In the case of GOC, lability is constant in the mixed layer
126	! It is computed only below the mixed layer depth
127	! ------------------------------------------------------------
128	!
129	IF( fsdept(ji,jj,jk) > zdep ) THEN
130	alphat = 0.
131	remint = 0.
132	!
133	zsizek1 = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
134	zsizek = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
135	!
136	IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
137	!
138	! The first level just below the mixed layer needs a
139	! specific treatment because lability is supposed constant
140	! everywhere within the mixed layer. This means that
141	! change in lability in the bottom part of the previous cell
142	! should not be computed
143	! ----------------------------------------------------------
144	!
145	! POC concentration is computed using the lagrangian
146	! framework. It is only used for the lability param
147	zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2 &
148	& * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
149	zpoc = max(0., zpoc)
150	!
151	DO jn = 1, jcpoc
152	!
153	! Lagrangian based algorithm. The fraction of each
154	! lability class is computed starting from the previous
155	! level
156	! -----------------------------------------------------
157	!
158	! the concentration of each lability class is calculated
159	! as the sum of the different sources and sinks
160	! Please note that production of new GOC experiences
161	! degradation
162	alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
163	& + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn) &
164	& * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2
165	alphat = alphat + alphag(ji,jj,jk,jn)
166	remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
167	END DO
168	ELSE
169	!
170	! standard algorithm in the rest of the water column
171	! See the comments in the previous block.
172	! ---------------------------------------------------
173	!
174	zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2 &
175	& * fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk) &
176	& * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
177	zpoc = max(0., zpoc)
178	!
179	DO jn = 1, jcpoc
180	alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek &
181	& + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1. &
182	& - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
183	& / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn)
184	alphat = alphat + alphag(ji,jj,jk,jn)
185	remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
186	END DO
187	ENDIF
188	DO jn = 1, jcpoc
189	! The contribution of each lability class at the current
190	! level is computed
191	alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn)
192	END DO
193	! Computation of the mean remineralisation rate
194	zremigoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint / ( alphat + rtrn) ))
195	ENDIF
196	ENDIF
197	END DO
198	END DO
199	END DO
200
201	DO jk = 1, jpkm1
202	DO jj = 1, jpj
203	DO ji = 1, jpi
204	zstep = xstep
205	# if defined key_degrad
206	zstep = zstep * facvol(ji,jj,jk)
207	# endif
208	! POC disaggregation by turbulence and bacterial activity.
209	! --------------------------------------------------------
210	zremig = zremigoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
211	zopoc2 = zremig * trb(ji,jj,jk,jpgoc)
212	orem(ji,jj,jk) = zopoc2
213	zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
214	zopon2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpgon)
215	zopop2 = xremipp / xremipc * zremig * trb(ji,jj,jk,jpgop)
216	zofer2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpbfe)
217
218	! Update the appropriate tracers trends
219	! -------------------------------------
220	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem3(ji,jj,jk)
221	tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + solgoc * zopon2
222	tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + solgoc * zopop2
223	tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + solgoc * zofer2
224	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc2
225	tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon2
226	tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop2
227	tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer2
228	zfolimi(ji,jj,jk) = zofer2
229	tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zopoc2 - zorem3(ji,jj,jk)
230	tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) - zopon2 * (1. + solgoc)
231	tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) - zopop2 * (1. + solgoc)
232	tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2 * (1. + solgoc)
233	END DO
234	END DO
235	END DO
236
237	IF(ln_ctl) THEN ! print mean trends (used for debugging)
238	WRITE(charout, FMT="('poc1')")
239	CALL prt_ctl_trc_info(charout)
240	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
241	ENDIF
242
243	! ------------------------------------------------------------------
244	! Lability parameterization for the small OM particles. This param
245	! is based on the same theoretical background as the big particles.
246	! However, because of its low sinking speed, lability is not supposed
247	! to be equal to its initial value (the value of the freshly produced
248	! organic matter). It is however uniform in the mixed layer.
249	! -------------------------------------------------------------------
250	!
251	totprod(:,:) = 0.
252	totthick(:,:) = 0.
253	totcons(:,:) = 0.
254	! intregrated production and consumption of POC in the mixed layer
255	! ----------------------------------------------------------------
256	!
257	DO jk = 1, jpkm1
258	DO jj = 1, jpj
259	DO ji = 1, jpi
260	zdep = hmld(ji,jj)
261	IF (tmask(ji,jj,jk) == 1. .AND. fsdept(ji,jj,jk) <= zdep ) THEN
262	totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
263	! The temperature effect is included here
264	totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
265	totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2 &
266	& / ( trb(ji,jj,jk,jppoc) + rtrn )
267	ENDIF
268	END DO
269	END DO
270	END DO
271
272	! Computation of the lability spectrum in the mixed layer. In the mixed
273	! layer, this spectrum is supposed to be uniform.
274	! ---------------------------------------------------------------------
275	DO jk = 1, jpkm1
276	DO jj = 1, jpj
277	DO ji = 1, jpi
278	IF (tmask(ji,jj,jk) == 1.) THEN
279	zdep = hmld(ji,jj)
280	alphat = 0.0
281	remint = 0.0
282	IF( fsdept(ji,jj,jk) <= zdep ) THEN
283	DO jn = 1, jcpoc
284	! For each lability class, the system is supposed to be
285	! at equilibrium: Prod - Sink - w alphap = 0.
286	alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn) &
287	& * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
288	alphat = alphat + alphap(ji,jj,jk,jn)
289	END DO
290	DO jn = 1, jcpoc
291	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
292	remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
293	END DO
294	! Mean remineralization rate in the mixed layer
295	zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint ))
296	ENDIF
297	ENDIF
298	END DO
299	END DO
300	END DO
301	!
302	! -----------------------------------------------------------------------
303	! The lability parameterization is used here. The code is here
304	! almost identical to what is done for big particles. The only difference
305	! is that an additional source from GOC to POC is included. This means
306	! that since we need the lability spectrum of GOC, GOC spectrum
307	! should be determined before.
308	! -----------------------------------------------------------------------
309	!
310	DO jk = 2, jpkm1
311	DO jj = 1, jpj
312	DO ji = 1, jpi
313	IF (tmask(ji,jj,jk) == 1.) THEN
314	zdep = hmld(ji,jj)
315	IF( fsdept(ji,jj,jk) > zdep ) THEN
316	alphat = 0.
317	remint = 0.
318	!
319	! the scale factors are corrected with temperature
320	zsizek1 = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
321	zsizek = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
322	!
323	! Special treatment of the level just below the MXL
324	! See the comments in the GOC section
325	! ---------------------------------------------------
326	!
327	IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
328	!
329	! Computation of the POC concentration using the
330	! lagrangian algorithm
331	zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2 &
332	& * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
333	zpoc = max(0., zpoc)
334	!
335	DO jn = 1, jcpoc
336	! computation of the lability spectrum applying the
337	! different sources and sinks
338	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
339	& + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
340	& / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn) &
341	& * zsizek ) )
342	alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
343	alphat = alphat + alphap(ji,jj,jk,jn)
344	END DO
345	ELSE
346	!
347	! Lability parameterization for the interior of the ocean
348	! This is very similar to what is done in the previous
349	! block
350	! --------------------------------------------------------
351	!
352	zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2 &
353	& * fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) &
354	& * rday / rfact2 * fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
355	zpoc = max(0., zpoc)
356	!
357	DO jn = 1, jcpoc
358	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) &
359	& * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn) &
360	& + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn) &
361	& / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) &
362	& * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) &
363	& * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1. &
364	& - exp( -reminp(jn) * zsizek ) )
365	alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
366	alphat = alphat + alphap(ji,jj,jk,jn)
367	END DO
368	ENDIF
369	! Normalization of the lability spectrum so that the
370	! integral is equal to 1
371	DO jn = 1, jcpoc
372	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
373	remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
374	END DO
375	! Mean remineralization rate in the water column
376	zremipoc(ji,jj,jk) = MIN(xremipc, MAX(0., remint ))
377	ENDIF
378	ENDIF
379	END DO
380	END DO
381	END DO
382	#endif
383
384	DO jk = 1, jpkm1
385	DO jj = 1, jpj
386	DO ji = 1, jpi
387	zstep = xstep
388	# if defined key_degrad
389	zstep = zstep * facvol(ji,jj,jk)
390	# endif
391	! POC disaggregation by turbulence and bacterial activity.
392	! --------------------------------------------------------
393	zremip = zremipoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
394	zopoc = zremip * trb(ji,jj,jk,jppoc)
395	orem(ji,jj,jk) = orem(ji,jj,jk) + zopoc
396	zopon = xremipn / xremipc * zremip * trb(ji,jj,jk,jppon)
397	zopop = xremipp / xremipc * zremip * trb(ji,jj,jk,jppop)
398	zofer = xremipn / xremipc * zremip * trb(ji,jj,jk,jpsfe)
399	#if defined key_kriest
400	zopoc2 = xremipc * zremip * trb(ji,jj,jk,jpnum)
401	#endif
402
403	! Update the appropriate tracers trends
404	! -------------------------------------
405	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zopoc
406	tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zopon
407	tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zopop
408	tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
409	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc
410	tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon
411	tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop
412	tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer
413	#if defined key_kriest
414	tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zopoc2
415	#endif
416	END DO
417	END DO
418	END DO
419
420	IF( ln_diatrc .AND. lk_iomput .AND. knt == nrdttrc ) THEN
421	zrfact2 = 1.e3 * rfact2r
422	CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) ) ! Remineralisation rate
423	CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) ) ! Remineralisation rate
424	CALL iom_put( "REMINF" , zfolimi(:,:,:) * tmask(:,:,:) * 1.e+3 * rfact2r ) ! Remineralisation rate
425	ENDIF
426	!
427	IF(ln_ctl) THEN ! print mean trends (used for debugging)
428	WRITE(charout, FMT="('poc1')")
429	CALL prt_ctl_trc_info(charout)
430	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
431	ENDIF
432	!
433	CALL wrk_dealloc( jpi, jpj, jpk, zwork1 )
434	CALL wrk_dealloc( jpi, jpj, totprod, totthick, totcons )
435	CALL wrk_dealloc( jpi, jpj, jpk, zremipoc, zremigoc, zorem3, zfolimi )
436	DEALLOCATE( alphag )
437	!
438	IF( nn_timing == 1 ) CALL timing_stop('p5z_poc')
439	!
440	END SUBROUTINE p5z_poc
441
442
443	SUBROUTINE p5z_poc_init
444	!!----------------------------------------------------------------------
445	!! * ROUTINE p5z_poc_init *
446	!!
447	!! ** Purpose : Initialization of remineralization parameters
448	!!
449	!! ** Method : Read the nampisrem namelist and check the parameters
450	!! called at the first timestep
451	!!
452	!! ** input : Namelist nampisrem
453	!!
454	!!----------------------------------------------------------------------
455	INTEGER :: jn
456	REAL(wp) :: remindelta, reminup, remindown
457	INTEGER :: ifault
458
459	NAMELIST/nampispoc/ xremipc, xremipn, xremipp, jcpoc, rshape
460	INTEGER :: ios ! Local integer output status for namelist read
461
462	REWIND( numnatp_ref ) ! Namelist nampispoc in reference namelist : Pisces remineralization
463	READ ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
464	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampispoc in reference namelist', lwp )
465
466	REWIND( numnatp_cfg ) ! Namelist nampispoc in configuration namelist : Pisces remineralization
467	READ ( numnatp_cfg, nampispoc, IOSTAT = ios, ERR = 902 )
468	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampispoc in configuration namelist', lwp )
469	IF(lwm) WRITE ( numonp, nampispoc )
470
471	IF(lwp) THEN ! control print
472	WRITE(numout,*) ' '
473	WRITE(numout,*) ' Namelist parameters for remineralization, nampispoc'
474	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
475	WRITE(numout,*) ' remineralisation rate of POC xremipc =', xremipc
476	WRITE(numout,*) ' remineralisation rate of PON xremipn =', xremipn
477	WRITE(numout,*) ' remineralisation rate of POP xremipp =', xremipp
478	WRITE(numout,*) ' Number of lability classes for POC jcpoc =', jcpoc
479	WRITE(numout,*) ' Shape factor of the gamma distribution rshape =', rshape
480	ENDIF
481	!
482	!
483	! Discretization along the lability space
484	! ---------------------------------------
485	!
486	ALLOCATE( alphan(jcpoc), reminp(jcpoc) )
487	ALLOCATE( alphap(jpi,jpj,jpk,jcpoc) )
488	!
489	IF (jcpoc > 1) THEN
490	!
491	remindelta = log(4. * 1000. ) / float(jcpoc-1)
492	reminup = 1./ 400. * exp(remindelta)
493	!
494	! Discretization based on incomplete gamma functions
495	! As incomplete gamma functions are not available in standard
496	! fortran 95, they have been coded as functions in this module (gamain)
497	! ---------------------------------------------------------------------
498	!
499	alphan(1) = gamain(reminup, rshape, ifault)
500	reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremipc / alphan(1)
501	DO jn = 2, jcpoc-1
502	reminup = 1./ 400. * exp(float(jn) * remindelta)
503	remindown = 1. / 400. * exp(float(jn-1) * remindelta)
504	alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault)
505	reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault)
506	reminp(jn) = reminp(jn) * xremipc / alphan(jn)
507	END DO
508	remindown = 1. / 400. * exp(float(jcpoc-1) * remindelta)
509	alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault)
510	reminp(jcpoc) = 1.0 - gamain(remindown, rshape+1.0, ifault)
511	reminp(jcpoc) = reminp(jcpoc) * xremipc / alphan(jcpoc)
512	ELSE
513	alphan(jcpoc) = 1.
514	reminp(jcpoc) = xremipc
515	ENDIF
516
517	DO jn = 1, jcpoc
518	alphap(:,:,:,jn) = alphan(jn)
519	END DO
520
521	END SUBROUTINE p5z_poc_init
522
523	REAL FUNCTION alngam( xvalue, ifault )
524
525	!*****************************************************************************80
526	!
527	!! ALNGAM computes the logarithm of the gamma function.
528	!
529	! Modified:
530	!
531	! 13 January 2008
532	!
533	! Author:
534	!
535	! Allan Macleod
536	! FORTRAN90 version by John Burkardt
537	!
538	! Reference:
539	!
540	! Allan Macleod,
541	! Algorithm AS 245,
542	! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
543	! Applied Statistics,
544	! Volume 38, Number 2, 1989, pages 397-402.
545	!
546	! Parameters:
547	!
548	! Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
549	!
550	! Output, integer ( kind = 4 ) IFAULT, error flag.
551	! 0, no error occurred.
552	! 1, XVALUE is less than or equal to 0.
553	! 2, XVALUE is too big.
554	!
555	! Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
556	!
557	implicit none
558
559	real(wp), parameter :: alr2pi = 0.918938533204673E+00
560	integer:: ifault
561	real(wp), dimension ( 9 ) :: r1 = (/ &
562	-2.66685511495E+00, &
563	-24.4387534237E+00, &
564	-21.9698958928E+00, &
565	11.1667541262E+00, &
566	3.13060547623E+00, &
567	0.607771387771E+00, &
568	11.9400905721E+00, &
569	31.4690115749E+00, &
570	15.2346874070E+00 /)
571	real(wp), dimension ( 9 ) :: r2 = (/ &
572	-78.3359299449E+00, &
573	-142.046296688E+00, &
574	137.519416416E+00, &
575	78.6994924154E+00, &
576	4.16438922228E+00, &
577	47.0668766060E+00, &
578	313.399215894E+00, &
579	263.505074721E+00, &
580	43.3400022514E+00 /)
581	real(wp), dimension ( 9 ) :: r3 = (/ &
582	-2.12159572323E+05, &
583	2.30661510616E+05, &
584	2.74647644705E+04, &
585	-4.02621119975E+04, &
586	-2.29660729780E+03, &
587	-1.16328495004E+05, &
588	-1.46025937511E+05, &
589	-2.42357409629E+04, &
590	-5.70691009324E+02 /)
591	real(wp), dimension ( 5 ) :: r4 = (/ &
592	0.279195317918525E+00, &
593	0.4917317610505968E+00, &
594	0.0692910599291889E+00, &
595	3.350343815022304E+00, &
596	6.012459259764103E+00 /)
597	real (wp) :: x
598	real (wp) :: x1
599	real (wp) :: x2
600	real (wp), parameter :: xlge = 5.10E+05
601	real (wp), parameter :: xlgst = 1.0E+30
602	real (wp) :: xvalue
603	real (wp) :: y
604
605	x = xvalue
606	alngam = 0.0E+00
607	!
608	! Check the input.
609	!
610	if ( xlgst <= x ) then
611	ifault = 2
612	return
613	end if
614	if ( x <= 0.0E+00 ) then
615	ifault = 1
616	return
617	end if
618	ifault = 0
619	!
620	! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
621	!
622	if ( x < 1.5E+00 ) then
623
624	if ( x < 0.5E+00 ) then
625	alngam = - log ( x )
626	y = x + 1.0E+00
627	!
628	! Test whether X < machine epsilon.
629	!
630	if ( y == 1.0E+00 ) then
631	return
632	end if
633
634	else
635
636	alngam = 0.0E+00
637	y = x
638	x = ( x - 0.5E+00 ) - 0.5E+00
639
640	end if
641
642	alngam = alngam + x * (((( &
643	r1(5) * y &
644	+ r1(4) ) * y &
645	+ r1(3) ) * y &
646	+ r1(2) ) * y &
647	+ r1(1) ) / (((( &
648	y &
649	+ r1(9) ) * y &
650	+ r1(8) ) * y &
651	+ r1(7) ) * y &
652	+ r1(6) )
653
654	return
655
656	end if
657	!
658	! Calculation for 1.5 <= X < 4.0.
659	!
660	if ( x < 4.0E+00 ) then
661
662	y = ( x - 1.0E+00 ) - 1.0E+00
663
664	alngam = y * (((( &
665	r2(5) * x &
666	+ r2(4) ) * x &
667	+ r2(3) ) * x &
668	+ r2(2) ) * x &
669	+ r2(1) ) / (((( &
670	x &
671	+ r2(9) ) * x &
672	+ r2(8) ) * x &
673	+ r2(7) ) * x &
674	+ r2(6) )
675	!
676	! Calculation for 4.0 <= X < 12.0.
677	!
678	else if ( x < 12.0E+00 ) then
679
680	alngam = (((( &
681	r3(5) * x &
682	+ r3(4) ) * x &
683	+ r3(3) ) * x &
684	+ r3(2) ) * x &
685	+ r3(1) ) / (((( &
686	x &
687	+ r3(9) ) * x &
688	+ r3(8) ) * x &
689	+ r3(7) ) * x &
690	+ r3(6) )
691	!
692	! Calculation for 12.0 <= X.
693	!
694	else
695
696	y = log ( x )
697	alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
698
699	if ( x <= xlge ) then
700
701	x1 = 1.0E+00 / x
702	x2 = x1 * x1
703
704	alngam = alngam + x1 * ( ( &
705	r4(3) * &
706	x2 + r4(2) ) * &
707	x2 + r4(1) ) / ( ( &
708	x2 + r4(5) ) * &
709	x2 + r4(4) )
710
711	end if
712	end if
713
714	END FUNCTION alngam
715
716	REAL FUNCTION gamain( x, p, ifault )
717
718	!*****************************************************************************80
719	!
720	!! GAMAIN computes the incomplete gamma ratio.
721	!
722	! Discussion:
723	!
724	! A series expansion is used if P > X or X <= 1. Otherwise, a
725	! continued fraction approximation is used.
726	!
727	! Modified:
728	!
729	! 17 January 2008
730	!
731	! Author:
732	!
733	! G Bhattacharjee
734	! FORTRAN90 version by John Burkardt
735	!
736	! Reference:
737	!
738	! G Bhattacharjee,
739	! Algorithm AS 32:
740	! The Incomplete Gamma Integral,
741	! Applied Statistics,
742	! Volume 19, Number 3, 1970, pages 285-287.
743	!
744	! Parameters:
745	!
746	! Input, real ( kind = 8 ) X, P, the parameters of the incomplete
747	! gamma ratio. 0 <= X, and 0 < P.
748	!
749	! Output, integer ( kind = 4 ) IFAULT, error flag.
750	! 0, no errors.
751	! 1, P <= 0.
752	! 2, X < 0.
753	! 3, underflow.
754	! 4, error return from the Log Gamma routine.
755	!
756	! Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
757	! gamma ratio.
758	!
759	implicit none
760	real (wp) a
761	real (wp), parameter :: acu = 1.0E-08
762	real (wp) an
763	real (wp) arg
764	real (wp) b
765	real (wp) dif
766	real (wp) factor
767	real (wp) g
768	real (wp) gin
769	integer i
770	integer ifault
771	real (wp), parameter :: oflo = 1.0E+37
772	real (wp) p
773	real (wp) pn(6)
774	real (wp) rn
775	real (wp) term
776	real (wp), parameter :: uflo = 1.0E-37
777	real (wp) x
778	!
779	! Check the input.
780	!
781	if ( p <= 0.0E+00 ) then
782	ifault = 1
783	gamain = 0.0E+00
784	return
785	end if
786
787	if ( x < 0.0E+00 ) then
788	ifault = 2
789	gamain = 0.0E+00
790	return
791	end if
792
793	if ( x == 0.0E+00 ) then
794	ifault = 0
795	gamain = 0.0E+00
796	return
797	end if
798
799	g = alngam ( p, ifault )
800
801	if ( ifault /= 0 ) then
802	ifault = 4
803	gamain = 0.0E+00
804	return
805	end if
806
807	arg = p * log ( x ) - x - g
808	if ( arg < log ( uflo ) ) then
809	ifault = 3
810	gamain = 0.0E+00
811	return
812	end if
813
814	ifault = 0
815	factor = exp ( arg )
816	!
817	! Calculation by series expansion.
818	!
819	if ( x <= 1.0E+00 .or. x < p ) then
820
821	gin = 1.0E+00
822	term = 1.0E+00
823	rn = p
824
825	do
826
827	rn = rn + 1.0E+00
828	term = term * x / rn
829	gin = gin + term
830
831	if ( term <= acu ) then
832	exit
833	end if
834
835	end do
836
837	gamain = gin * factor / p
838	return
839
840	end if
841	!
842	! Calculation by continued fraction.
843	!
844	a = 1.0E+00 - p
845	b = a + x + 1.0E+00
846	term = 0.0E+00
847
848	pn(1) = 1.0E+00
849	pn(2) = x
850	pn(3) = x + 1.0E+00
851	pn(4) = x * b
852
853	gin = pn(3) / pn(4)
854	do
855
856	a = a + 1.0E+00
857	b = b + 2.0E+00
858	term = term + 1.0E+00
859	an = a * term
860	do i = 1, 2
861	pn(i+4) = b * pn(i+2) - an * pn(i)
862	end do
863
864	if ( pn(6) /= 0.0E+00 ) then
865
866	rn = pn(5) / pn(6)
867	dif = abs ( gin - rn )
868	!
869	! Absolute error tolerance satisfied?
870	!
871	if ( dif <= acu ) then
872	!
873	! Relative error tolerance satisfied?
874	!
875	if ( dif <= acu * rn ) then
876	gamain = 1.0E+00 - factor * gin
877	exit
878	end if
879
880	end if
881
882	gin = rn
883
884	end if
885
886	do i = 1, 4
887	pn(i) = pn(i+2)
888	end do
889	if ( oflo <= abs ( pn(5) ) ) then
890
891	do i = 1, 4
892	pn(i) = pn(i) / oflo
893	end do
894
895	end if
896
897	end do
898
899	END FUNCTION gamain
900
901	#else
902	!!======================================================================
903	!! Dummy module : No PISCES bio-model
904	!!======================================================================
905	CONTAINS
906	SUBROUTINE p5z_poc ! Empty routine
907	END SUBROUTINE p5z_poc
908	#endif
909
910	!!======================================================================
911	END MODULE p5zpoc

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

Download in other formats: