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.

p4zsink.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Context Navigation

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 6841

Last change on this file since 6841 was 6841, checked in by aumont, 8 years ago
Various bug fixes + explicit gamma function for lability
File size: 33.4 KB

Line
1	MODULE p4zsink
2	!!======================================================================
3	!! * MODULE p4zsink *
4	!! TOP : PISCES vertical flux of particulate matter due to gravitational sinking
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) Change aggregation formula
9	!! 3.5 ! 2012-07 (O. Aumont) Introduce potential time-splitting
10	!!----------------------------------------------------------------------
11	#if defined key_pisces
12	!!----------------------------------------------------------------------
13	!! p4z_sink : Compute vertical flux of particulate matter due to gravitational sinking
14	!! p4z_sink_init : Unitialisation of sinking speed parameters
15	!! p4z_sink_alloc : Allocate sinking speed variables
16	!!----------------------------------------------------------------------
17	USE oce_trc ! shared variables between ocean and passive tracers
18	USE trc ! passive tracers common variables
19	USE sms_pisces ! PISCES Source Minus Sink variables
20	USE prtctl_trc ! print control for debugging
21	USE iom ! I/O manager
22	USE lib_mpp
23
24	IMPLICIT NONE
25	PRIVATE
26
27	PUBLIC p4z_sink ! called in p4zbio.F90
28	PUBLIC p4z_sink_init ! called in trcsms_pisces.F90
29	PUBLIC p4z_sink_alloc
30
31	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsbio3 !: POC sinking speed
32	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsbio4 !: GOC sinking speed
33	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wscal !: Calcite and BSi sinking speeds
34
35	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinking, sinking2 !: POC sinking fluxes
36	! ! (different meanings depending on the parameterization)
37	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkcal, sinksil !: CaCO3 and BSi sinking fluxes
38	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer !: Small BFe sinking fluxes
39	#if ! defined key_kriest
40	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfer2 !: Big iron sinking fluxes
41	#endif
42	#if defined key_ligand
43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sinkfep !: Fep sinking fluxes
44	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wsfep
45	#endif
46
47	INTEGER :: ik100
48
49	#if defined key_kriest
50	REAL(wp) :: xkr_sfact !: Sinking factor
51	REAL(wp) :: xkr_stick !: Stickiness
52	REAL(wp) :: xkr_nnano !: Nbr of cell in nano size class
53	REAL(wp) :: xkr_ndiat !: Nbr of cell in diatoms size class
54	REAL(wp) :: xkr_nmicro !: Nbr of cell in microzoo size class
55	REAL(wp) :: xkr_nmeso !: Nbr of cell in mesozoo size class
56	REAL(wp) :: xkr_naggr !: Nbr of cell in aggregates size class
57
58	REAL(wp) :: xkr_frac
59
60	REAL(wp), PUBLIC :: xkr_dnano !: Size of particles in nano pool
61	REAL(wp), PUBLIC :: xkr_ddiat !: Size of particles in diatoms pool
62	REAL(wp), PUBLIC :: xkr_dmicro !: Size of particles in microzoo pool
63	REAL(wp), PUBLIC :: xkr_dmeso !: Size of particles in mesozoo pool
64	REAL(wp), PUBLIC :: xkr_daggr !: Size of particles in aggregates pool
65	REAL(wp), PUBLIC :: xkr_wsbio_min !: min vertical particle speed
66	REAL(wp), PUBLIC :: xkr_wsbio_max !: max vertical particle speed
67
68	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: xnumm !: maximum number of particles in aggregates
69	#endif
70
71	!!* Substitution
72	# include "top_substitute.h90"
73	!!----------------------------------------------------------------------
74	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
75	!! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
76	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
77	!!----------------------------------------------------------------------
78	CONTAINS
79
80	#if ! defined key_kriest
81	!!----------------------------------------------------------------------
82	!! 'standard sinking parameterisation' ???
83	!!----------------------------------------------------------------------
84
85	SUBROUTINE p4z_sink ( kt, knt )
86	!!---------------------------------------------------------------------
87	!! * ROUTINE p4z_sink *
88	!!
89	!! ** Purpose : Compute vertical flux of particulate matter due to
90	!! gravitational sinking
91	!!
92	!! ** Method : - ???
93	!!---------------------------------------------------------------------
94	INTEGER, INTENT(in) :: kt, knt
95	INTEGER :: ji, jj, jk, jit
96	INTEGER :: iiter1, iiter2
97	REAL(wp) :: zfact, zwsmax, zmax
98	CHARACTER (len=25) :: charout
99	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
100	REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d
101	!!---------------------------------------------------------------------
102	!
103	IF( nn_timing == 1 ) CALL timing_start('p4z_sink')
104
105	! Initialization of some global variables
106	! ---------------------------------------
107	prodpoc(:,:,:) = 0.
108	conspoc(:,:,:) = 0.
109	prodgoc(:,:,:) = 0.
110	consgoc(:,:,:) = 0.
111	!
112	! Sinking speeds of detritus is increased with depth as shown
113	! by data and from the coagulation theory
114	! -----------------------------------------------------------
115	DO jk = 1, jpkm1
116	DO jj = 1, jpj
117	DO ji = 1,jpi
118	zmax = MAX( heup_01(ji,jj), hmld(ji,jj) )
119	zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale
120	wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
121	END DO
122	END DO
123	END DO
124
125	! limit the values of the sinking speeds to avoid numerical instabilities
126	wsbio3(:,:,:) = wsbio
127	wscal(:,:,:) = wsbio4(:,:,:)
128	#if defined key_ligand
129	wsfep (:,:,:) = wfep
130	#endif
131	!
132	! OA This is (I hope) a temporary solution for the problem that may
133	! OA arise in specific situation where the CFL criterion is broken
134	! OA for vertical sedimentation of particles. To avoid this, a time
135	! OA splitting algorithm has been coded. A specific maximum
136	! OA iteration number is provided and may be specified in the namelist
137	! OA This is to avoid very large iteration number when explicit free
138	! OA surface is used (for instance). When niter?max is set to 1,
139	! OA this computation is skipped. The crude old threshold method is
140	! OA then applied. This also happens when niter exceeds nitermax.
141	IF( MAX( niter1max, niter2max ) == 1 ) THEN
142	iiter1 = 1
143	iiter2 = 1
144	ELSE
145	iiter1 = 1
146	iiter2 = 1
147	DO jk = 1, jpkm1
148	DO jj = 1, jpj
149	DO ji = 1, jpi
150	IF( tmask(ji,jj,jk) == 1) THEN
151	zwsmax = 0.5 * fse3t(ji,jj,jk) / xstep
152	iiter1 = MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
153	iiter2 = MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
154	ENDIF
155	END DO
156	END DO
157	END DO
158	IF( lk_mpp ) THEN
159	CALL mpp_max( iiter1 )
160	CALL mpp_max( iiter2 )
161	ENDIF
162	iiter1 = MIN( iiter1, niter1max )
163	iiter2 = MIN( iiter2, niter2max )
164	ENDIF
165
166	DO jk = 1,jpkm1
167	DO jj = 1, jpj
168	DO ji = 1, jpi
169	IF( tmask(ji,jj,jk) == 1 ) THEN
170	zwsmax = 0.5 * fse3t(ji,jj,jk) / xstep
171	wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
172	wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
173	#if defined key_ligand
174	wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
175	#endif
176	ENDIF
177	END DO
178	END DO
179	END DO
180
181	! Initializa to zero all the sinking arrays
182	! -----------------------------------------
183	sinking (:,:,:) = 0.e0
184	sinking2(:,:,:) = 0.e0
185	sinkcal (:,:,:) = 0.e0
186	sinkfer (:,:,:) = 0.e0
187	sinksil (:,:,:) = 0.e0
188	sinkfer2(:,:,:) = 0.e0
189	#if defined key_ligand
190	sinkfep(:,:,:) = 0.e0
191	#endif
192
193	! Compute the sedimentation term using p4zsink2 for all the sinking particles
194	! -----------------------------------------------------
195	DO jit = 1, iiter1
196	CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
197	CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
198	#if defined key_ligand
199	CALL p4z_sink2( wsfep , sinkfep , jpfep, iiter1 )
200	#endif
201	END DO
202
203	DO jit = 1, iiter2
204	CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
205	CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
206	CALL p4z_sink2( wscal, sinksil , jpgsi, iiter2 )
207	CALL p4z_sink2( wscal, sinkcal , jpcal, iiter2 )
208	END DO
209
210	! Total carbon export per year
211	IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &
212	& t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
213	!
214	IF( lk_iomput ) THEN
215	IF( knt == nrdttrc ) THEN
216	CALL wrk_alloc( jpi, jpj, zw2d )
217	CALL wrk_alloc( jpi, jpj, jpk, zw3d )
218	zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
219	!
220	IF( iom_use( "EPC100" ) ) THEN
221	zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
222	CALL iom_put( "EPC100" , zw2d )
223	ENDIF
224	IF( iom_use( "EPFE100" ) ) THEN
225	zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
226	CALL iom_put( "EPFE100" , zw2d )
227	ENDIF
228	IF( iom_use( "EPCAL100" ) ) THEN
229	zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
230	CALL iom_put( "EPCAL100" , zw2d )
231	ENDIF
232	IF( iom_use( "EPSI100" ) ) THEN
233	zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
234	CALL iom_put( "EPSI100" , zw2d )
235	ENDIF
236	IF( iom_use( "EXPC" ) ) THEN
237	zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
238	CALL iom_put( "EXPC" , zw3d )
239	ENDIF
240	IF( iom_use( "EXPFE" ) ) THEN
241	zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron
242	CALL iom_put( "EXPFE" , zw3d )
243	ENDIF
244	IF( iom_use( "EXPCAL" ) ) THEN
245	zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
246	CALL iom_put( "EXPCAL" , zw3d )
247	ENDIF
248	IF( iom_use( "EXPSI" ) ) THEN
249	zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
250	CALL iom_put( "EXPSI" , zw3d )
251	ENDIF
252	IF( iom_use( "tcexp" ) ) CALL iom_put( "tcexp" , t_oce_co2_exp * zfact ) ! molC/s
253	!
254	CALL wrk_dealloc( jpi, jpj, zw2d )
255	CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
256	ENDIF
257	ELSE
258	IF( ln_diatrc ) THEN
259	zfact = 1.e3 * rfact2r
260	trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
261	trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
262	trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
263	trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)
264	trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
265	trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
266	ENDIF
267	ENDIF
268	!
269	IF(ln_ctl) THEN ! print mean trends (used for debugging)
270	WRITE(charout, FMT="('sink')")
271	CALL prt_ctl_trc_info(charout)
272	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
273	ENDIF
274	!
275	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink')
276	!
277	END SUBROUTINE p4z_sink
278
279	SUBROUTINE p4z_sink_init
280	!!----------------------------------------------------------------------
281	!! * ROUTINE p4z_sink_init *
282	!!----------------------------------------------------------------------
283	INTEGER :: jk
284
285	ik100 = 10 ! last level where depth less than 100 m
286	DO jk = jpkm1, 1, -1
287	IF( gdept_1d(jk) > 100. ) ik100 = jk - 1
288	END DO
289	IF (lwp) WRITE(numout,*)
290	IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1
291	IF (lwp) WRITE(numout,*)
292	!
293	t_oce_co2_exp = 0._wp
294	!
295	END SUBROUTINE p4z_sink_init
296
297	#else
298	!!----------------------------------------------------------------------
299	!! 'Kriest sinking parameterisation' key_kriest ???
300	!!----------------------------------------------------------------------
301
302	SUBROUTINE p4z_sink ( kt, knt )
303	!!---------------------------------------------------------------------
304	!! * ROUTINE p4z_sink *
305	!!
306	!! ** Purpose : Compute vertical flux of particulate matter due to
307	!! gravitational sinking - Kriest parameterization
308	!!
309	!! ** Method : - ???
310	!!---------------------------------------------------------------------
311	!
312	INTEGER, INTENT(in) :: kt, knt
313	!
314	INTEGER :: ji, jj, jk, jit, niter1, niter2
315	REAL(wp) :: znum , zeps, zfm, zgm, zsm
316	REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
317	REAL(wp) :: zval1, zval2, zval3
318	REAL(wp) :: zfact
319	INTEGER :: ik1
320	CHARACTER (len=25) :: charout
321	REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d
322	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
323	REAL(wp), POINTER, DIMENSION(:,: ) :: zw2d
324	!!---------------------------------------------------------------------
325	!
326	IF( nn_timing == 1 ) CALL timing_start('p4z_sink')
327	!
328	CALL wrk_alloc( jpi, jpj, jpk, znum3d )
329	!
330	! Initialisation of variables used to compute Sinking Speed
331	! ---------------------------------------------------------
332
333	znum3d(:,:,:) = 0.e0
334	zval1 = 1. + xkr_zeta
335	zval2 = 1. + xkr_zeta + xkr_eta
336	zval3 = 1. + xkr_eta
337
338	! Computation of the vertical sinking speed : Kriest et Evans, 2000
339	! -----------------------------------------------------------------
340
341	DO jk = 1, jpkm1
342	DO jj = 1, jpj
343	DO ji = 1, jpi
344	IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
345	znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
346	! -------------- To avoid sinking speed over 50 m/day -------
347	znum = MIN( xnumm(jk), znum )
348	znum = MAX( 1.1 , znum )
349	znum3d(ji,jj,jk) = znum
350	!------------------------------------------------------------
351	zeps = ( zval1 * znum - 1. )/ ( znum - 1. )
352	zfm = xkr_frac**( 1. - zeps )
353	zgm = xkr_frac**( zval1 - zeps )
354	zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
355	zdiv1 = zeps - zval3
356	wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &
357	& - xkr_wsbio_max * zgm * xkr_eta / zdiv
358	wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &
359	& - xkr_wsbio_max * zfm * xkr_eta / zdiv1
360	IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
361	ENDIF
362	END DO
363	END DO
364	END DO
365
366	wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )
367	#if defined key_ligand
368	wsfep (:,:,:) = wfep
369	#endif
370
371	! INITIALIZE TO ZERO ALL THE SINKING ARRAYS
372	! -----------------------------------------
373
374	sinking (:,:,:) = 0.e0
375	sinking2(:,:,:) = 0.e0
376	sinkcal (:,:,:) = 0.e0
377	sinkfer (:,:,:) = 0.e0
378	sinksil (:,:,:) = 0.e0
379	#if defined key_ligand
380	sinkfep(:,:,:) = 0.e0
381	#endif
382
383	! Compute the sedimentation term using p4zsink2 for all the sinking particles
384	! -----------------------------------------------------
385
386	niter1 = niter1max
387	niter2 = niter2max
388
389	DO jit = 1, niter1
390	CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )
391	CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )
392	CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
393	CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
394	#if defined key_ligand
395	CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 )
396	#endif
397	END DO
398
399	DO jit = 1, niter2
400	CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )
401	END DO
402
403	IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc ) ) &
404	& t_oce_co2_exp = glob_sum( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )
405	!
406	IF( lk_iomput ) THEN
407	IF( knt == nrdttrc ) THEN
408	CALL wrk_alloc( jpi, jpj, zw2d )
409	CALL wrk_alloc( jpi, jpj, jpk, zw3d )
410	zfact = 1.e+3 * rfact2r ! conversion from mol/l/kt to mol/m3/s
411	!
412	IF( iom_use( "EPC100" ) ) THEN
413	zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m
414	CALL iom_put( "EPC100" , zw2d )
415	ENDIF
416	IF( iom_use( "EPN100" ) ) THEN
417	zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ?
418	CALL iom_put( "EPN100" , zw2d )
419	ENDIF
420	IF( iom_use( "EPCAL100" ) ) THEN
421	zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
422	CALL iom_put( "EPCAL100" , zw2d )
423	ENDIF
424	IF( iom_use( "EPSI100" ) ) THEN
425	zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
426	CALL iom_put( "EPSI100" , zw2d )
427	ENDIF
428	IF( iom_use( "EXPC" ) ) THEN
429	zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
430	CALL iom_put( "EXPC" , zw3d )
431	ENDIF
432	IF( iom_use( "EXPN" ) ) THEN
433	zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
434	CALL iom_put( "EXPN" , zw3d )
435	ENDIF
436	IF( iom_use( "EXPCAL" ) ) THEN
437	zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
438	CALL iom_put( "EXPCAL" , zw3d )
439	ENDIF
440	IF( iom_use( "EXPSI" ) ) THEN
441	zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
442	CALL iom_put( "EXPSI" , zw3d )
443	ENDIF
444	IF( iom_use( "XNUM" ) ) THEN
445	zw3d(:,:,:) = znum3d(:,:,:) * tmask(:,:,:) ! Number of particles on aggregats
446	CALL iom_put( "XNUM" , zw3d )
447	ENDIF
448	IF( iom_use( "WSC" ) ) THEN
449	zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles
450	CALL iom_put( "WSC" , zw3d )
451	ENDIF
452	IF( iom_use( "WSN" ) ) THEN
453	zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number
454	CALL iom_put( "WSN" , zw3d )
455	ENDIF
456	!
457	CALL wrk_dealloc( jpi, jpj, zw2d )
458	CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
459	ELSE
460	IF( ln_diatrc ) THEN
461	zfact = 1.e3 * rfact2r
462	trc2d(:,: ,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
463	trc2d(:,: ,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
464	trc2d(:,: ,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
465	trc2d(:,: ,jp_pcs0_2d + 7) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
466	trc2d(:,: ,jp_pcs0_2d + 8) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
467	trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:) * zfact * tmask(:,:,:)
468	trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:) * zfact * tmask(:,:,:)
469	trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:) * zfact * tmask(:,:,:)
470	trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:) * zfact * tmask(:,:,:)
471	trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d (:,:,:) * tmask(:,:,:)
472	trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3 (:,:,:) * tmask(:,:,:)
473	trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4 (:,:,:) * tmask(:,:,:)
474	ENDIF
475	ENDIF
476
477	!
478	IF(ln_ctl) THEN ! print mean trends (used for debugging)
479	WRITE(charout, FMT="('sink')")
480	CALL prt_ctl_trc_info(charout)
481	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
482	ENDIF
483	!
484	CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
485	!
486	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink')
487	!
488	END SUBROUTINE p4z_sink
489
490
491	SUBROUTINE p4z_sink_init
492	!!----------------------------------------------------------------------
493	!! * ROUTINE p4z_sink_init *
494	!!
495	!! ** Purpose : Initialization of sinking parameters
496	!! Kriest parameterization only
497	!!
498	!! ** Method : Read the nampiskrs namelist and check the parameters
499	!! called at the first timestep
500	!!
501	!! ** input : Namelist nampiskrs
502	!!----------------------------------------------------------------------
503	INTEGER :: jk, jn, kiter
504	INTEGER :: ios ! Local integer output status for namelist read
505	REAL(wp) :: znum, zdiv
506	REAL(wp) :: zws, zwr, zwl,wmax, znummax
507	REAL(wp) :: zmin, zmax, zl, zr, xacc
508	!
509	NAMELIST/nampiskrs/ xkr_sfact, xkr_stick , &
510	& xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr
511	!!----------------------------------------------------------------------
512	!
513	IF( nn_timing == 1 ) CALL timing_start('p4z_sink_init')
514	!
515
516	REWIND( numnatp_ref ) ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest
517	READ ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)
518	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )
519
520	REWIND( numnatp_cfg ) ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest
521	READ ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )
522	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )
523	IF(lwm) WRITE ( numonp, nampiskrs )
524
525	IF(lwp) THEN
526	WRITE(numout,*)
527	WRITE(numout,*) ' Namelist : nampiskrs'
528	WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact
529	WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick
530	WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano
531	WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat
532	WRITE(numout,*) ' Nbr of cell in microzoo size class xkr_nmicro = ', xkr_nmicro
533	WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso
534	WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr
535	ENDIF
536
537
538	! max and min vertical particle speed
539	xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
540	xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
541	IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
542
543	!
544	! effect of the sizes of the different living pools on particle numbers
545	! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
546	! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
547	! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
548	! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
549	! doc aggregates = 1um
550	! ----------------------------------------------------------
551
552	xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
553	xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
554	xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )
555	xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
556	xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
557
558	!!---------------------------------------------------------------------
559	!! 'key_kriest' ???
560	!!---------------------------------------------------------------------
561	! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
562	! Search of the maximum number of particles in aggregates for each k-level.
563	! Bissection Method
564	!--------------------------------------------------------------------
565	IF (lwp) THEN
566	WRITE(numout,*)
567	WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'
568	ENDIF
569
570	xacc = 0.001_wp
571	kiter = 50
572	zmin = 1.10_wp
573	zmax = xkr_mass_max / xkr_mass_min
574	xkr_frac = zmax
575
576	DO jk = 1,jpk
577	zl = zmin
578	zr = zmax
579	wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2
580	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
581	znum = zl - 1.
582	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
583	& - ( xkr_wsbio_max * xkr_eta * znum * &
584	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
585	& - wmax
586
587	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
588	znum = zr - 1.
589	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
590	& - ( xkr_wsbio_max * xkr_eta * znum * &
591	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
592	& - wmax
593	iflag: DO jn = 1, kiter
594	IF ( zwl == 0._wp ) THEN ; znummax = zl
595	ELSEIF( zwr == 0._wp ) THEN ; znummax = zr
596	ELSE
597	znummax = ( zr + zl ) / 2.
598	zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
599	znum = znummax - 1.
600	zws = xkr_wsbio_min * xkr_zeta / zdiv &
601	& - ( xkr_wsbio_max * xkr_eta * znum * &
602	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
603	& - wmax
604	IF( zws * zwl < 0. ) THEN ; zr = znummax
605	ELSE ; zl = znummax
606	ENDIF
607	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
608	znum = zl - 1.
609	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
610	& - ( xkr_wsbio_max * xkr_eta * znum * &
611	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
612	& - wmax
613
614	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
615	znum = zr - 1.
616	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
617	& - ( xkr_wsbio_max * xkr_eta * znum * &
618	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
619	& - wmax
620	!
621	IF ( ABS ( zws ) <= xacc ) EXIT iflag
622	!
623	ENDIF
624	!
625	END DO iflag
626
627	xnumm(jk) = znummax
628	IF (lwp) WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
629	!
630	END DO
631	!
632	ik100 = 10 ! last level where depth less than 100 m
633	DO jk = jpkm1, 1, -1
634	IF( gdept_1d(jk) > 100. ) ik100 = jk - 1
635	END DO
636	IF (lwp) WRITE(numout,*)
637	IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ', ik100 + 1
638	IF (lwp) WRITE(numout,*)
639	!
640	t_oce_co2_exp = 0._wp
641	!
642	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink_init')
643	!
644	END SUBROUTINE p4z_sink_init
645
646	#endif
647
648	SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
649	!!---------------------------------------------------------------------
650	!! * ROUTINE p4z_sink2 *
651	!!
652	!! ** Purpose : Compute the sedimentation terms for the various sinking
653	!! particles. The scheme used to compute the trends is based
654	!! on MUSCL.
655	!!
656	!! ** Method : - this ROUTINE compute not exactly the advection but the
657	!! transport term, i.e. div(u*tra).
658	!!---------------------------------------------------------------------
659	!
660	INTEGER , INTENT(in ) :: jp_tra ! tracer index index
661	INTEGER , INTENT(in ) :: kiter ! number of iterations for time-splitting
662	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed
663	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe
664	!!
665	INTEGER :: ji, jj, jk, jn
666	REAL(wp) :: zigma,zew,zign, zflx, zstep
667	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb
668	!!---------------------------------------------------------------------
669	!
670	IF( nn_timing == 1 ) CALL timing_start('p4z_sink2')
671	!
672	! Allocate temporary workspace
673	CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
674
675	zstep = rfact2 / FLOAT( kiter ) / 2.
676
677	ztraz(:,:,:) = 0.e0
678	zakz (:,:,:) = 0.e0
679	ztrb (:,:,:) = trb(:,:,:,jp_tra)
680
681	DO jk = 1, jpkm1
682	zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)
683	END DO
684	zwsink2(:,:,1) = 0.e0
685	IF( lk_degrad ) THEN
686	zwsink2(:,:,:) = zwsink2(:,:,:) * facvol(:,:,:)
687	ENDIF
688
689
690	! Vertical advective flux
691	DO jn = 1, 2
692	! first guess of the slopes interior values
693	DO jk = 2, jpkm1
694	ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
695	END DO
696	ztraz(:,:,1 ) = 0.0
697	ztraz(:,:,jpk) = 0.0
698
699	! slopes
700	DO jk = 2, jpkm1
701	DO jj = 1,jpj
702	DO ji = 1, jpi
703	zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
704	zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
705	END DO
706	END DO
707	END DO
708
709	! Slopes limitation
710	DO jk = 2, jpkm1
711	DO jj = 1, jpj
712	DO ji = 1, jpi
713	zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * &
714	& MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
715	END DO
716	END DO
717	END DO
718
719	! vertical advective flux
720	DO jk = 1, jpkm1
721	DO jj = 1, jpj
722	DO ji = 1, jpi
723	zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
724	zew = zwsink2(ji,jj,jk+1)
725	psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
726	END DO
727	END DO
728	END DO
729	!
730	! Boundary conditions
731	psinkflx(:,:,1 ) = 0.e0
732	psinkflx(:,:,jpk) = 0.e0
733
734	DO jk=1,jpkm1
735	DO jj = 1,jpj
736	DO ji = 1, jpi
737	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
738	trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
739	END DO
740	END DO
741	END DO
742
743	ENDDO
744
745	DO jk = 1,jpkm1
746	DO jj = 1,jpj
747	DO ji = 1, jpi
748	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
749	ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
750	END DO
751	END DO
752	END DO
753
754	trb(:,:,:,jp_tra) = ztrb(:,:,:)
755	psinkflx(:,:,:) = 2. * psinkflx(:,:,:)
756	!
757	CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
758	!
759	IF( nn_timing == 1 ) CALL timing_stop('p4z_sink2')
760	!
761	END SUBROUTINE p4z_sink2
762
763
764	INTEGER FUNCTION p4z_sink_alloc()
765	!!----------------------------------------------------------------------
766	!! * ROUTINE p4z_sink_alloc *
767	!!----------------------------------------------------------------------
768	ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4 (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) , &
769	& sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk) , &
770	& sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk) , &
771	#if defined key_kriest
772	& xnumm(jpk) , &
773	#else
774	& sinkfer2(jpi,jpj,jpk) , &
775	#endif
776	#if defined key_ligand
777	& wsfep(jpi,jpj,jpk) , sinkfep(jpi,jpj,jpk) , &
778	#endif
779	& sinkfer(jpi,jpj,jpk) , STAT=p4z_sink_alloc )
780	!
781	IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
782	!
783	END FUNCTION p4z_sink_alloc
784
785	#else
786	!!======================================================================
787	!! Dummy module : No PISCES bio-model
788	!!======================================================================
789	CONTAINS
790	SUBROUTINE p4z_sink ! Empty routine
791	END SUBROUTINE p4z_sink
792	#endif
793
794	!!======================================================================
795	END MODULE p4zsink

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

Download in other formats: