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 NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/TOP/PISCES/P4Z – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 13176

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

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

Download in other formats: