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

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

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

Download in other formats: