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.

p4zpoc.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90 @ 9125

Last change on this file since 9125 was 9125, checked in by timgraham, 6 years ago
Removed wrk_arrays from whole code. No change in SETTE results from this.
File size: 33.5 KB

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

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

Download in other formats: