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/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Context Navigation

source: branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90 @ 7192

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

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

Download in other formats: