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

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

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

Download in other formats: