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 @ 6453

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

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

Download in other formats: