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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/PISCES/P4Z – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 14219

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

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

Download in other formats: