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

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

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

Download in other formats: