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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
Property svn:keywords set to `Id`
File size: 34.7 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	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	!!----------------------------------------------------------------------
40	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
41	!! $Id$
42	!! Software governed by the CeCILL license (see ./LICENSE)
43	!!----------------------------------------------------------------------
44	CONTAINS
45
46	SUBROUTINE p4z_poc( kt, knt, Kbb, Kmm, Krhs )
47	!!---------------------------------------------------------------------
48	!! * ROUTINE p4z_poc *
49	!!
50	!! ** Purpose : Compute remineralization of organic particles
51	!!
52	!! ** Method : - ???
53	!!---------------------------------------------------------------------
54	INTEGER, INTENT(in) :: kt, knt ! ocean time step and ???
55	INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices
56	!
57	INTEGER :: ji, jj, jk, jn
58	REAL(wp) :: zremip, zremig, zdep, zorem, zorem2, zofer
59	REAL(wp) :: zopon, zopop, zopoc, zopoc2, zopon2, zopop2
60	REAL(wp) :: zsizek, zsizek1, alphat, remint, solgoc, zpoc
61	REAL(wp) :: zofer2, zofer3
62	REAL(wp) :: zrfact2
63	CHARACTER (len=25) :: charout
64	REAL(wp), DIMENSION(jpi,jpj ) :: totprod, totthick, totcons
65	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
66	REAL(wp), DIMENSION(jpi,jpj,jpk,jcpoc) :: alphag
67	!!---------------------------------------------------------------------
68	!
69	IF( ln_timing ) CALL timing_start('p4z_poc')
70	!
71	! Initialization of local variables
72	! ---------------------------------
73
74	! Here we compute the GOC -> POC rate due to the shrinking
75	! of the fecal pellets/aggregates as a result of bacterial
76	! solubilization
77	! This is based on a fractal dimension of 2.56 and a spectral
78	! slope of -3.6 (identical to what is used in p4zsink to compute
79	! aggregation
80	solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
81
82	! Initialisation of temprary arrys
83	IF( ln_p4z ) THEN
84	zremipoc(:,:,:) = xremip
85	zremigoc(:,:,:) = xremip
86	ELSE ! ln_p5z
87	zremipoc(:,:,:) = xremipc
88	zremigoc(:,:,:) = xremipc
89	ENDIF
90	zorem3(:,:,:) = 0.
91	orem (:,:,:) = 0.
92	ztremint(:,:,:) = 0.
93	zfolimi (:,:,:) = 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	ztremint(:,:,:) = zremigoc(:,:,:)
106	DO jk = 2, jpkm1
107	DO jj = 1, jpj
108	DO ji = 1, jpi
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 DO
187	END DO
188	END DO
189
190	IF( ln_p4z ) THEN ; zremigoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
191	ELSE ; zremigoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
192	ENDIF
193
194	IF( ln_p4z ) THEN
195	DO jk = 1, jpkm1
196	DO jj = 1, jpj
197	DO ji = 1, jpi
198	! POC disaggregation by turbulence and bacterial activity.
199	! --------------------------------------------------------
200	zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
201	zorem2 = zremig * tr(ji,jj,jk,jpgoc,Kbb)
202	orem(ji,jj,jk) = zorem2
203	zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
204	zofer2 = zremig * tr(ji,jj,jk,jpbfe,Kbb)
205	zofer3 = zremig * solgoc * tr(ji,jj,jk,jpbfe,Kbb)
206
207	! -------------------------------------
208	tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
209	tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zorem2 - zorem3(ji,jj,jk)
210	tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zofer3
211	tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 - zofer3
212	tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem2
213	tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
214	zfolimi(ji,jj,jk) = zofer2
215	END DO
216	END DO
217	END DO
218	ELSE
219	DO jk = 1, jpkm1
220	DO jj = 1, jpj
221	DO ji = 1, jpi
222	! POC disaggregation by turbulence and bacterial activity.
223	! --------------------------------------------------------
224	zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
225	zopoc2 = zremig * tr(ji,jj,jk,jpgoc,Kbb)
226	orem(ji,jj,jk) = zopoc2
227	zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
228	zopon2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpgon,Kbb)
229	zopop2 = xremipp / xremipc * zremig * tr(ji,jj,jk,jpgop,Kbb)
230	zofer2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpbfe,Kbb)
231
232	! -------------------------------------
233	tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
234	tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + solgoc * zopon2
235	tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + solgoc * zopop2
236	tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + solgoc * zofer2
237	tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc2
238	tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon2
239	tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop2
240	tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
241	tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zopoc2 - zorem3(ji,jj,jk)
242	tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) - zopon2 * (1. + solgoc)
243	tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) - zopop2 * (1. + solgoc)
244	tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 * (1. + solgoc)
245	zfolimi(ji,jj,jk) = zofer2
246	END DO
247	END DO
248	END DO
249	ENDIF
250
251	IF(ln_ctl) THEN ! print mean trends (used for debugging)
252	WRITE(charout, FMT="('poc1')")
253	CALL prt_ctl_trc_info(charout)
254	CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
255	ENDIF
256
257	! ------------------------------------------------------------------
258	! Lability parameterization for the small OM particles. This param
259	! is based on the same theoretical background as the big particles.
260	! However, because of its low sinking speed, lability is not supposed
261	! to be equal to its initial value (the value of the freshly produced
262	! organic matter). It is however uniform in the mixed layer.
263	! -------------------------------------------------------------------
264	!
265	totprod (:,:) = 0.
266	totthick(:,:) = 0.
267	totcons (:,:) = 0.
268	! intregrated production and consumption of POC in the mixed layer
269	! ----------------------------------------------------------------
270	!
271	DO jk = 1, jpkm1
272	DO jj = 1, jpj
273	DO ji = 1, jpi
274	zdep = hmld(ji,jj)
275	IF (tmask(ji,jj,jk) == 1. .AND. gdept(ji,jj,jk,Kmm) <= zdep ) THEN
276	totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2
277	! The temperature effect is included here
278	totthick(ji,jj) = totthick(ji,jj) + e3t(ji,jj,jk,Kmm)* tgfunc(ji,jj,jk)
279	totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2 &
280	& / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn )
281	ENDIF
282	END DO
283	END DO
284	END DO
285
286	! Computation of the lability spectrum in the mixed layer. In the mixed
287	! layer, this spectrum is supposed to be uniform.
288	! ---------------------------------------------------------------------
289	ztremint(:,:,:) = zremipoc(:,:,:)
290	DO jk = 1, jpkm1
291	DO jj = 1, jpj
292	DO ji = 1, jpi
293	IF (tmask(ji,jj,jk) == 1.) THEN
294	zdep = hmld(ji,jj)
295	alphat = 0.0
296	remint = 0.0
297	IF( gdept(ji,jj,jk,Kmm) <= zdep ) THEN
298	DO jn = 1, jcpoc
299	! For each lability class, the system is supposed to be
300	! at equilibrium: Prod - Sink - w alphap = 0.
301	alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn) &
302	& * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
303	alphat = alphat + alphap(ji,jj,jk,jn)
304	END DO
305	DO jn = 1, jcpoc
306	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
307	remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
308	END DO
309	! Mean remineralization rate in the mixed layer
310	ztremint(ji,jj,jk) = MAX( 0., remint )
311	ENDIF
312	ENDIF
313	END DO
314	END DO
315	END DO
316	!
317	IF( ln_p4z ) THEN ; zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
318	ELSE ; zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
319	ENDIF
320
321	! -----------------------------------------------------------------------
322	! The lability parameterization is used here. The code is here
323	! almost identical to what is done for big particles. The only difference
324	! is that an additional source from GOC to POC is included. This means
325	! that since we need the lability spectrum of GOC, GOC spectrum
326	! should be determined before.
327	! -----------------------------------------------------------------------
328	!
329	DO jk = 2, jpkm1
330	DO jj = 1, jpj
331	DO ji = 1, jpi
332	IF (tmask(ji,jj,jk) == 1.) THEN
333	zdep = hmld(ji,jj)
334	IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
335	alphat = 0.
336	remint = 0.
337	!
338	! the scale factors are corrected with temperature
339	zsizek1 = e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
340	zsizek = e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
341	!
342	! Special treatment of the level just below the MXL
343	! See the comments in the GOC section
344	! ---------------------------------------------------
345	!
346	IF ( gdept(ji,jj,jk-1,Kmm) <= zdep ) THEN
347	!
348	! Computation of the POC concentration using the
349	! lagrangian algorithm
350	zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk) * rday / rfact2 &
351	& * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
352	zpoc = max(0., zpoc)
353	!
354	DO jn = 1, jcpoc
355	! computation of the lability spectrum applying the
356	! different sources and sinks
357	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
358	& + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
359	& / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn) &
360	& * zsizek ) )
361	alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
362	alphat = alphat + alphap(ji,jj,jk,jn)
363	END DO
364	ELSE
365	!
366	! Lability parameterization for the interior of the ocean
367	! This is very similar to what is done in the previous
368	! block
369	! --------------------------------------------------------
370	!
371	zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk-1) * rday / rfact2 &
372	& * e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk) &
373	& * rday / rfact2 * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
374	zpoc = max(0., zpoc)
375	!
376	DO jn = 1, jcpoc
377	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) &
378	& * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn) &
379	& + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn) &
380	& / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) &
381	& * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) &
382	& * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1. &
383	& - exp( -reminp(jn) * zsizek ) )
384	alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
385	alphat = alphat + alphap(ji,jj,jk,jn)
386	END DO
387	ENDIF
388	! Normalization of the lability spectrum so that the
389	! integral is equal to 1
390	DO jn = 1, jcpoc
391	alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
392	remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
393	END DO
394	! Mean remineralization rate in the water column
395	ztremint(ji,jj,jk) = MAX( 0., remint )
396	ENDIF
397	ENDIF
398	END DO
399	END DO
400	END DO
401
402	IF( ln_p4z ) THEN ; zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
403	ELSE ; zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
404	ENDIF
405
406	IF( ln_p4z ) THEN
407	DO jk = 1, jpkm1
408	DO jj = 1, jpj
409	DO ji = 1, jpi
410	IF (tmask(ji,jj,jk) == 1.) THEN
411	! POC disaggregation by turbulence and bacterial activity.
412	! --------------------------------------------------------
413	zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
414	zorem = zremip * tr(ji,jj,jk,jppoc,Kbb)
415	zofer = zremip * tr(ji,jj,jk,jpsfe,Kbb)
416
417	tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem
418	orem(ji,jj,jk) = orem(ji,jj,jk) + zorem
419	tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer
420	tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zorem
421	tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
422	zfolimi(ji,jj,jk) = zfolimi(ji,jj,jk) + 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 * tr(ji,jj,jk,jppoc,Kbb)
435	orem(ji,jj,jk) = orem(ji,jj,jk) + zopoc
436	zopon = xremipn / xremipc * zremip * tr(ji,jj,jk,jppon,Kbb)
437	zopop = xremipp / xremipc * zremip * tr(ji,jj,jk,jppop,Kbb)
438	zofer = xremipn / xremipc * zremip * tr(ji,jj,jk,jpsfe,Kbb)
439
440	tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zopoc
441	tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zopon
442	tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) - zopop
443	tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
444	tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc
445	tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon
446	tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop
447	tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer
448	zfolimi(ji,jj,jk) = zfolimi(ji,jj,jk) + 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	CALL iom_put( "REMINF" , zfolimi(:,:,:) * tmask(:,:,:) * 1.e+9 * zrfact2 ) ! Remineralisation rate
460	ENDIF
461	ENDIF
462
463	IF(ln_ctl) THEN ! print mean trends (used for debugging)
464	WRITE(charout, FMT="('poc2')")
465	CALL prt_ctl_trc_info(charout)
466	CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
467	ENDIF
468	!
469	!
470	IF( ln_timing ) CALL timing_stop('p4z_poc')
471	!
472	END SUBROUTINE p4z_poc
473
474
475	SUBROUTINE p4z_poc_init
476	!!----------------------------------------------------------------------
477	!! * ROUTINE p4z_poc_init *
478	!!
479	!! ** Purpose : Initialization of remineralization parameters
480	!!
481	!! ** Method : Read the nampispoc namelist and check the parameters
482	!! called at the first timestep
483	!!
484	!! ** input : Namelist nampispoc
485	!!----------------------------------------------------------------------
486	INTEGER :: jn ! dummy loop index
487	INTEGER :: ios, ifault ! Local integer
488	REAL(wp):: remindelta, reminup, remindown
489	!!
490	NAMELIST/nampispoc/ xremip , jcpoc , rshape, &
491	& xremipc, xremipn, xremipp
492	!!----------------------------------------------------------------------
493	!
494	IF(lwp) THEN
495	WRITE(numout,*)
496	WRITE(numout,*) 'p4z_poc_init : Initialization of remineralization parameters'
497	WRITE(numout,*) '~~~~~~~~~~~~'
498	ENDIF
499	!
500	REWIND( numnatp_ref ) ! Namelist nampisrem in reference namelist : Pisces remineralization
501	READ ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
502	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampispoc in reference namelist' )
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' )
506	IF(lwm) WRITE( numonp, nampispoc )
507
508	IF(lwp) THEN ! control print
509	WRITE(numout,*) ' Namelist : nampispoc'
510	IF( ln_p4z ) THEN
511	WRITE(numout,*) ' remineralisation rate of POC xremip =', xremip
512	ELSE
513	WRITE(numout,*) ' remineralisation rate of POC xremipc =', xremipc
514	WRITE(numout,*) ' remineralisation rate of PON xremipn =', xremipn
515	WRITE(numout,*) ' remineralisation rate of POP xremipp =', xremipp
516	ENDIF
517	WRITE(numout,*) ' Number of lability classes for POC jcpoc =', jcpoc
518	WRITE(numout,*) ' Shape factor of the gamma distribution rshape =', rshape
519	ENDIF
520	!
521	! Discretization along the lability space
522	! ---------------------------------------
523	!
524	ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , 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
562	REAL FUNCTION alngam( xvalue, ifault )
563	!*****************************************************************************80
564	!
565	!! ALNGAM computes the logarithm of the gamma function.
566	!
567	! Modified: 13 January 2008
568	!
569	! Author : Allan Macleod
570	! FORTRAN90 version by John Burkardt
571	!
572	! Reference:
573	! Allan Macleod, Algorithm AS 245,
574	! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
575	! Applied Statistics,
576	! Volume 38, Number 2, 1989, pages 397-402.
577	!
578	! Parameters:
579	!
580	! Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
581	!
582	! Output, integer ( kind = 4 ) IFAULT, error flag.
583	! 0, no error occurred.
584	! 1, XVALUE is less than or equal to 0.
585	! 2, XVALUE is too big.
586	!
587	! Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
588	!*****************************************************************************80
589	implicit none
590
591	real(wp), parameter :: alr2pi = 0.918938533204673E+00
592	integer:: ifault
593	real(wp), dimension ( 9 ) :: r1 = (/ &
594	-2.66685511495E+00, &
595	-24.4387534237E+00, &
596	-21.9698958928E+00, &
597	11.1667541262E+00, &
598	3.13060547623E+00, &
599	0.607771387771E+00, &
600	11.9400905721E+00, &
601	31.4690115749E+00, &
602	15.2346874070E+00 /)
603	real(wp), dimension ( 9 ) :: r2 = (/ &
604	-78.3359299449E+00, &
605	-142.046296688E+00, &
606	137.519416416E+00, &
607	78.6994924154E+00, &
608	4.16438922228E+00, &
609	47.0668766060E+00, &
610	313.399215894E+00, &
611	263.505074721E+00, &
612	43.3400022514E+00 /)
613	real(wp), dimension ( 9 ) :: r3 = (/ &
614	-2.12159572323E+05, &
615	2.30661510616E+05, &
616	2.74647644705E+04, &
617	-4.02621119975E+04, &
618	-2.29660729780E+03, &
619	-1.16328495004E+05, &
620	-1.46025937511E+05, &
621	-2.42357409629E+04, &
622	-5.70691009324E+02 /)
623	real(wp), dimension ( 5 ) :: r4 = (/ &
624	0.279195317918525E+00, &
625	0.4917317610505968E+00, &
626	0.0692910599291889E+00, &
627	3.350343815022304E+00, &
628	6.012459259764103E+00 /)
629	real (wp) :: x
630	real (wp) :: x1
631	real (wp) :: x2
632	real (wp), parameter :: xlge = 5.10E+05
633	real (wp), parameter :: xlgst = 1.0E+30
634	real (wp) :: xvalue
635	real (wp) :: y
636
637	x = xvalue
638	alngam = 0.0E+00
639	!
640	! Check the input.
641	!
642	if ( xlgst <= x ) then
643	ifault = 2
644	return
645	end if
646	if ( x <= 0.0E+00 ) then
647	ifault = 1
648	return
649	end if
650
651	ifault = 0
652	!
653	! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
654	!
655	if ( x < 1.5E+00 ) then
656
657	if ( x < 0.5E+00 ) then
658	alngam = - log ( x )
659	y = x + 1.0E+00
660	!
661	! Test whether X < machine epsilon.
662	!
663	if ( y == 1.0E+00 ) then
664	return
665	end if
666
667	else
668
669	alngam = 0.0E+00
670	y = x
671	x = ( x - 0.5E+00 ) - 0.5E+00
672
673	end if
674
675	alngam = alngam + x * (((( &
676	r1(5) * y &
677	+ r1(4) ) * y &
678	+ r1(3) ) * y &
679	+ r1(2) ) * y &
680	+ r1(1) ) / (((( &
681	y &
682	+ r1(9) ) * y &
683	+ r1(8) ) * y &
684	+ r1(7) ) * y &
685	+ r1(6) )
686
687	return
688
689	end if
690	!
691	! Calculation for 1.5 <= X < 4.0.
692	!
693	if ( x < 4.0E+00 ) then
694
695	y = ( x - 1.0E+00 ) - 1.0E+00
696
697	alngam = y * (((( &
698	r2(5) * x &
699	+ r2(4) ) * x &
700	+ r2(3) ) * x &
701	+ r2(2) ) * x &
702	+ r2(1) ) / (((( &
703	x &
704	+ r2(9) ) * x &
705	+ r2(8) ) * x &
706	+ r2(7) ) * x &
707	+ r2(6) )
708	!
709	! Calculation for 4.0 <= X < 12.0.
710	!
711	else if ( x < 12.0E+00 ) then
712
713	alngam = (((( &
714	r3(5) * x &
715	+ r3(4) ) * x &
716	+ r3(3) ) * x &
717	+ r3(2) ) * x &
718	+ r3(1) ) / (((( &
719	x &
720	+ r3(9) ) * x &
721	+ r3(8) ) * x &
722	+ r3(7) ) * x &
723	+ r3(6) )
724	!
725	! Calculation for 12.0 <= X.
726	!
727	else
728
729	y = log ( x )
730	alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
731
732	if ( x <= xlge ) then
733
734	x1 = 1.0E+00 / x
735	x2 = x1 * x1
736
737	alngam = alngam + x1 * ( ( &
738	r4(3) * &
739	x2 + r4(2) ) * &
740	x2 + r4(1) ) / ( ( &
741	x2 + r4(5) ) * &
742	x2 + r4(4) )
743
744	end if
745
746	end if
747
748	END FUNCTION alngam
749
750
751	REAL FUNCTION gamain( x, p, ifault )
752	!*****************************************************************************80
753	!
754	!! GAMAIN computes the incomplete gamma ratio.
755	!
756	! Discussion:
757	!
758	! A series expansion is used if P > X or X <= 1. Otherwise, a
759	! continued fraction approximation is used.
760	!
761	! Modified:
762	!
763	! 17 January 2008
764	!
765	! Author:
766	!
767	! G Bhattacharjee
768	! FORTRAN90 version by John Burkardt
769	!
770	! Reference:
771	!
772	! G Bhattacharjee,
773	! Algorithm AS 32:
774	! The Incomplete Gamma Integral,
775	! Applied Statistics,
776	! Volume 19, Number 3, 1970, pages 285-287.
777	!
778	! Parameters:
779	!
780	! Input, real ( kind = 8 ) X, P, the parameters of the incomplete
781	! gamma ratio. 0 <= X, and 0 < P.
782	!
783	! Output, integer ( kind = 4 ) IFAULT, error flag.
784	! 0, no errors.
785	! 1, P <= 0.
786	! 2, X < 0.
787	! 3, underflow.
788	! 4, error return from the Log Gamma routine.
789	!
790	! Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
791	! gamma ratio.
792	!
793	implicit none
794
795	real (wp) a
796	real (wp), parameter :: acu = 1.0E-08
797	real (wp) an
798	real (wp) arg
799	real (wp) b
800	real (wp) dif
801	real (wp) factor
802	real (wp) g
803	real (wp) gin
804	integer i
805	integer ifault
806	real (wp), parameter :: oflo = 1.0E+37
807	real (wp) p
808	real (wp) pn(6)
809	real (wp) rn
810	real (wp) term
811	real (wp), parameter :: uflo = 1.0E-37
812	real (wp) x
813	!
814	! Check the input.
815	!
816	if ( p <= 0.0E+00 ) then
817	ifault = 1
818	gamain = 0.0E+00
819	return
820	end if
821
822	if ( x < 0.0E+00 ) then
823	ifault = 2
824	gamain = 0.0E+00
825	return
826	end if
827
828	if ( x == 0.0E+00 ) then
829	ifault = 0
830	gamain = 0.0E+00
831	return
832	end if
833
834	g = alngam ( p, ifault )
835
836	if ( ifault /= 0 ) then
837	ifault = 4
838	gamain = 0.0E+00
839	return
840	end if
841
842	arg = p * log ( x ) - x - g
843
844	if ( arg < log ( uflo ) ) then
845	ifault = 3
846	gamain = 0.0E+00
847	return
848	end if
849
850	ifault = 0
851	factor = exp ( arg )
852	!
853	! Calculation by series expansion.
854	!
855	if ( x <= 1.0E+00 .or. x < p ) then
856
857	gin = 1.0E+00
858	term = 1.0E+00
859	rn = p
860
861	do
862
863	rn = rn + 1.0E+00
864	term = term * x / rn
865	gin = gin + term
866
867	if ( term <= acu ) then
868	exit
869	end if
870
871	end do
872
873	gamain = gin * factor / p
874	return
875
876	end if
877	!
878	! Calculation by continued fraction.
879	!
880	a = 1.0E+00 - p
881	b = a + x + 1.0E+00
882	term = 0.0E+00
883
884	pn(1) = 1.0E+00
885	pn(2) = x
886	pn(3) = x + 1.0E+00
887	pn(4) = x * b
888
889	gin = pn(3) / pn(4)
890
891	do
892
893	a = a + 1.0E+00
894	b = b + 2.0E+00
895	term = term + 1.0E+00
896	an = a * term
897	do i = 1, 2
898	pn(i+4) = b * pn(i+2) - an * pn(i)
899	end do
900
901	if ( pn(6) /= 0.0E+00 ) then
902
903	rn = pn(5) / pn(6)
904	dif = abs ( gin - rn )
905	!
906	! Absolute error tolerance satisfied?
907	!
908	if ( dif <= acu ) then
909	!
910	! Relative error tolerance satisfied?
911	!
912	if ( dif <= acu * rn ) then
913	gamain = 1.0E+00 - factor * gin
914	exit
915	end if
916
917	end if
918
919	gin = rn
920
921	end if
922
923	do i = 1, 4
924	pn(i) = pn(i+2)
925	end do
926	if ( oflo <= abs ( pn(5) ) ) then
927
928	do i = 1, 4
929	pn(i) = pn(i) / oflo
930	end do
931
932	end if
933
934	end do
935
936	END FUNCTION gamain
937
938	!!======================================================================
939	END MODULE p4zpoc

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

Download in other formats: