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 trunk/NEMO/TOP_SRC/PISCES – NEMO

Context Navigation

source: trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90 @ 1075

Last change on this file since 1075 was 1073, checked in by cetlod, 16 years ago
update PISCES model, see ticket:190
Property svn:executable set to ``*
File size: 28.0 KB

Line
1	MODULE p4zsink
2	!!======================================================================
3	!! * MODULE p4zsink *
4	!! TOP : PISCES Compute 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	#if defined key_pisces
9	!!----------------------------------------------------------------------
10	!! p4z_sink : Compute vertical flux of particulate matter due to gravitational sinking
11	!!----------------------------------------------------------------------
12	USE trc
13	USE oce_trc !
14	USE trp_trc
15	USE sms_pisces
16	USE prtctl_trc
17
18
19	IMPLICIT NONE
20	PRIVATE
21
22	PUBLIC p4z_sink ! called in p4zbio.F90
23
24	!! * Shared module variables
25	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !:
26	wsbio3, wsbio4, & !: POC and GOC sinking speeds
27	wscal !: Calcite and BSi sinking speeds
28
29	!! * Module variables
30	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & !:
31	sinking, sinking2, & !: POC sinking fluxes (different meanings depending on the parameterization
32	sinkcal, sinksil, & !: CaCO3 and BSi sinking fluxes
33	sinkfer !: Small BFe sinking flux
34
35	#if defined key_kriest
36	REAL(wp) :: &
37	xkr_sfact = 250. , & !: Sinking factor
38	xkr_stick = 0.2 , & !: Stickiness
39	xkr_nnano = 2.337 , & !: Nbr of cell in nano size class
40	xkr_ndiat = 3.718 , & !: Nbr of cell in diatoms size class
41	xkr_nmeso = 7.147 , & !: Nbr of cell in mesozoo size class
42	xkr_naggr = 9.877 !: Nbr of cell in aggregates size class
43
44	REAL(wp) :: &
45	xkr_frac
46
47	REAL(wp), PUBLIC :: &
48	xkr_dnano , & !: Size of particles in nano pool
49	xkr_ddiat , & !: Size of particles in diatoms pool
50	xkr_dmeso , & !: Size of particles in mesozoo pool
51	xkr_daggr , & !: Size of particles in aggregates pool
52	xkr_wsbio_min , & !: min vertical particle speed
53	xkr_wsbio_max !: max vertical particle speed
54
55	REAL(wp), PUBLIC, DIMENSION(jpk) :: & !:
56	xnumm !: maximum number of particles in aggregates
57
58	#endif
59
60	#if ! defined key_kriest
61	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & !:
62	sinkfer2 !: Big Fe sinking flux
63	#endif
64
65	!!* Substitution
66	# include "domzgr_substitute.h90"
67	!!----------------------------------------------------------------------
68	!! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
69	!! $Header:$
70	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
71	!!----------------------------------------------------------------------
72
73	CONTAINS
74
75	#if defined key_kriest
76
77	SUBROUTINE p4z_sink ( kt, jnt )
78	!!---------------------------------------------------------------------
79	!! * ROUTINE p4z_sink *
80	!!
81	!! ** Purpose : Compute vertical flux of particulate matter due to
82	!! gravitational sinking - Kriest parameterization
83	!!
84	!! ** Method : - ???
85	!!---------------------------------------------------------------------
86
87	INTEGER, INTENT(in) :: kt, jnt
88	INTEGER :: ji, jj, jk
89	INTEGER :: iksed
90	REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
91	REAL(wp) :: zagg , zaggdoc, znumdoc
92	REAL(wp) :: znum , zeps, zfm, zgm, zsm
93	REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
94	REAL(wp) :: zval1, zval2, zval3, zval4
95	REAL(wp) :: zstep
96	#if defined key_trc_dia3d
97	REAL(wp) :: zrfact2
98	#endif
99	REAL(wp), DIMENSION(jpi,jpj,jpk) :: znum3d
100	CHARACTER (len=25) :: charout
101
102	!!---------------------------------------------------------------------
103
104	IF( ( kt * jnt ) == nittrc000 ) CALL p4z_sink_init ! Initialization (first time-step only)
105
106	zstep = rfact2 / rjjss ! Time step duration for biology
107
108
109	! Initialisation of variables used to compute Sinking Speed
110	! ---------------------------------------------------------
111
112	znum3d(:,:,:) = 0.e0
113	iksed = 10
114	zval1 = 1. + xkr_zeta
115	zval2 = 1. + xkr_zeta + xkr_eta
116	zval3 = 1. + xkr_eta
117
118	! Computation of the vertical sinking speed : Kriest et Evans, 2000
119	! -----------------------------------------------------------------
120
121	DO jk = 1, jpkm1
122	DO jj = 1, jpj
123	DO ji = 1, jpi
124	IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
125	znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
126	! -------------- To avoid sinking speed over 50 m/day -------
127	znum = MIN( xnumm(jk), znum )
128	znum = MAX( 1.1 , znum )
129	znum3d(ji,jj,jk) = znum
130	!------------------------------------------------------------
131	zeps = ( zval1 * znum - 1. )/ ( znum - 1. )
132	zfm = xkr_frac**( 1. - zeps )
133	zgm = xkr_frac**( zval1 - zeps )
134	zdiv = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
135	zdiv1 = zeps - zval3
136	wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv &
137	& - xkr_wsbio_max * zgm * xkr_eta / zdiv
138	wsbio4(ji,jj,jk) = xkr_wsbio_min * ( zeps-1. ) / zdiv1 &
139	& - xkr_wsbio_max * zfm * xkr_eta / zdiv1
140	IF( znum == 1.1) wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
141	ENDIF
142	END DO
143	END DO
144	END DO
145
146	wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
147
148
149	! INITIALIZE TO ZERO ALL THE SINKING ARRAYS
150	! -----------------------------------------
151
152	sinking (:,:,:) = 0.e0
153	sinking2(:,:,:) = 0.e0
154	sinkcal (:,:,:) = 0.e0
155	sinkfer (:,:,:) = 0.e0
156	sinksil (:,:,:) = 0.e0
157
158	! Compute the sedimentation term using p4zsink2 for all
159	! the sinking particles
160	! -----------------------------------------------------
161
162	CALL p4z_sink2( wsbio3, sinking , jppoc )
163	CALL p4z_sink2( wsbio4, sinking2, jpnum )
164	CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
165	CALL p4z_sink2( wscal , sinksil , jpdsi )
166	CALL p4z_sink2( wscal , sinkcal , jpcal )
167
168	! Exchange between organic matter compartments due to
169	! coagulation/disaggregation
170	! ---------------------------------------------------
171
172	zval1 = 1. + xkr_zeta
173	zval2 = 1. + xkr_eta
174	zval3 = 3. + xkr_eta
175	zval4 = 4. + xkr_eta
176
177	DO jk = 1,jpkm1
178	DO jj = 1,jpj
179	DO ji = 1,jpi
180	IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
181
182	znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
183	! -------------- To avoid sinking speed over 50 m/day -------
184	znum = min(xnumm(jk),znum)
185	znum = MAX( 1.1,znum)
186	!------------------------------------------------------------
187	zeps = ( zval1 * znum - 1.) / ( znum - 1.)
188	zdiv = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
189	zdiv1 = MAX( 1.e-4, ABS( zeps - 4. ) ) * SIGN( 1., zeps - 4. )
190	zdiv2 = zeps - 2.
191	zdiv3 = zeps - 3.
192	zdiv4 = zeps - zval2
193	zdiv5 = 2.* zeps - zval4
194	zfm = xkr_frac**( 1.- zeps )
195	zsm = xkr_frac**xkr_eta
196
197	! Part I : Coagulation dependant on turbulence
198	! ----------------------------------------------
199
200	zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2 &
201	& * 2.( (zfm-1.)(zfmxkr_mass_max3-xkr_mass_min*3) &
202	& * (zeps-1)/zdiv1 + 3.(zfmxkr_mass_max-xkr_mass_min) &
203	& * (zfmxkr_mass_max2-xkr_mass_min*2) &
204	& * (zeps-1.)*2/(zdiv2zdiv3)) &
205	# if defined key_off_degrad
206	& *facvol(ji,jj,jk) &
207	# endif
208	& )
209
210	zagg2 = ( 20.163trn(ji,jj,jk,jpnum)*2zfm* &
211	& ((xkr_mass_max*3+3.(xkr_mass_max**2 &
212	& xkr_mass_min(zeps-1.)/zdiv2 &
213	& +xkr_mass_maxxkr_mass_min2(zeps-1.)/zdiv3) &
214	& +xkr_mass_min*3(zeps-1)/zdiv1) &
215	& -zfmxkr_mass_max3(1.+3.*((zeps-1.)/ &
216	& (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1)) &
217	# if defined key_off_degrad
218	& *facvol(ji,jj,jk) &
219	# endif
220	& )
221
222	zagg3 = ( 0.163trn(ji,jj,jk,jpnum)2zfm*28. * xkr_mass_max**3 &
223	# if defined key_off_degrad
224	& *facvol(ji,jj,jk) &
225	# endif
226	& )
227
228	zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
229
230	! Aggregation of small into large particles
231	! Part II : Differential settling
232	! ----------------------------------------------
233
234	zagg4 = ( 2.3.1410.125trn(ji,jj,jk,jpnum)2 &
235	& xkr_wsbio_min(zeps-1.)*2 &
236	& (xkr_mass_min2((1.-zsmzfm)/(zdiv3zdiv4) &
237	& -(1.-zfm)/(zdiv*(zeps-1.)))- &
238	& ((zfmzfmxkr_mass_max*2zsm-xkr_mass_min**2) &
239	& xkr_eta)/(zdivzdiv3*zdiv5) ) &
240	# if defined key_off_degrad
241	& *facvol(ji,jj,jk) &
242	# endif
243	& )
244
245	zagg5 = ( 2.3.1410.125trn(ji,jj,jk,jpnum)*2 &
246	& (zeps-1.)zfm*xkr_wsbio_min &
247	& (zsm(xkr_mass_min*2-zfmxkr_mass_max**2) &
248	& /zdiv3-(xkr_mass_min*2-zfmzsmxkr_mass_max*2) &
249	& /zdiv) &
250	# if defined key_off_degrad
251	& *facvol(ji,jj,jk) &
252	# endif
253	& )
254
255	zaggsi = ( zagg4 + zagg5 ) * zstep / 10.
256
257	zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
258
259	! Aggregation of DOC to small particles
260	! --------------------------------------
261
262	zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc) &
263	& + 1018. * trn(ji,jj,jk,jppoc) ) * zstep &
264	# if defined key_off_degrad
265	& * facvol(ji,jj,jk) &
266	# endif
267	& * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
268
269	znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
270	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc
271	tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zaggdoc * znumdoc - zagg
272	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc
273
274	ENDIF
275	END DO
276	END DO
277	END DO
278
279	# if defined key_trc_dia3d
280	zrfact2 = 1.e3 * rfact2r
281	trc2d(:,:, 5) = sinking (:,:,iksed+1) * zrfact2
282	trc2d(:,:, 6) = sinking2(:,:,iksed+1) * zrfact2
283	trc2d(:,:, 7) = sinkfer (:,:,iksed+1) * zrfact2
284	trc2d(:,:, 9) = sinksil (:,:,iksed+1) * zrfact2
285	trc2d(:,:,10) = sinkcal (:,:,iksed+1) * zrfact2
286	trc3d(:,:,:,12) = sinking (:,:,:) * zrfact2
287	trc3d(:,:,:,13) = sinking2(:,:,:) * zrfact2
288	trc3d(:,:,:,14) = sinksil (:,:,:) * zrfact2
289	trc3d(:,:,:,15) = sinkcal (:,:,:) * zrfact2
290	trc3d(:,:,:,16) = znum3d (:,:,:)
291	trc3d(:,:,:,17) = wsbio3 (:,:,:)
292	trc3d(:,:,:,18) = wsbio4 (:,:,:)
293	# endif
294	!
295	IF(ln_ctl) THEN ! print mean trends (used for debugging)
296	WRITE(charout, FMT="('sink')")
297	CALL prt_ctl_trc_info(charout)
298	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
299	ENDIF
300
301	END SUBROUTINE p4z_sink
302
303	SUBROUTINE p4z_sink_init
304	!!----------------------------------------------------------------------
305	!! * ROUTINE p4z_sink_init *
306	!!
307	!! ** Purpose : Initialization of sinking parameters
308	!! Kriest parameterization only
309	!!
310	!! ** Method : Read the natkriest namelist and check the parameters
311	!! called at the first timestep (nittrc000)
312	!!
313	!! ** input : Namelist natkriest
314	!!
315	!!----------------------------------------------------------------------
316	INTEGER :: jk, jn, kiter
317	REAL(wp) :: znum, zdiv
318	REAL(wp) :: zws, zwr, zwl,wmax, znummax
319	REAL(wp) :: zmin, zmax, zl, zr, xacc
320
321	NAMELIST/natkrsize/ xkr_sfact, xkr_stick , &
322	& xkr_nnano, xkr_ndiat, xkr_nmeso, xkr_naggr
323
324	!!----------------------------------------------------------------------
325	! ! natkriest : kriest parameters
326	! ! -----------------------------
327	REWIND( numnat ) ! read natkriest
328	READ ( numnat, natkrsize )
329
330	IF(lwp) THEN
331	WRITE(numout,*)
332	WRITE(numout,*) ' Namelist : natkrsize'
333	WRITE(numout,*) ' Sinking factor xkr_sfact = ', xkr_sfact
334	WRITE(numout,*) ' Stickiness xkr_stick = ', xkr_stick
335	WRITE(numout,*) ' Nbr of cell in nano size class xkr_nnano = ', xkr_nnano
336	WRITE(numout,*) ' Nbr of cell in diatoms size class xkr_ndiat = ', xkr_ndiat
337	WRITE(numout,*) ' Nbr of cell in mesozoo size class xkr_nmeso = ', xkr_nmeso
338	WRITE(numout,*) ' Nbr of cell in aggregates size class xkr_naggr = ', xkr_naggr
339	ENDIF
340
341
342	! max and min vertical particle speed
343	xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
344	xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
345	WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
346
347	!
348	! effect of the sizes of the different living pools on particle numbers
349	! nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
350	! diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
351	! mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
352	! aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
353	! doc aggregates = 1um
354	! ----------------------------------------------------------
355
356	xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
357	xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
358	xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
359	xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
360
361	!!---------------------------------------------------------------------
362	!! 'key_kriest' ???
363	!!---------------------------------------------------------------------
364	! COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
365	! Search of the maximum number of particles in aggregates for each k-level.
366	! Bissection Method
367	!--------------------------------------------------------------------
368	WRITE(numout,*)
369	WRITE(numout,*)' kriest : Compute maximum number of particles in aggregates'
370
371	xacc = 0.001
372	kiter = 50
373	zmin = 1.10
374	zmax = xkr_mass_max / xkr_mass_min
375	xkr_frac = zmax
376
377	DO jk = 1,jpk
378	zl = zmin
379	zr = zmax
380	wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
381	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
382	znum = zl - 1.
383	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
384	& - ( xkr_wsbio_max * xkr_eta * znum * &
385	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
386	& - wmax
387
388	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
389	znum = zr - 1.
390	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
391	& - ( xkr_wsbio_max * xkr_eta * znum * &
392	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
393	& - wmax
394	iflag: DO jn = 1, kiter
395	IF( zwl == 0.e0 ) THEN
396	znummax = zl
397	ELSE IF ( zwr == 0.e0 ) THEN
398	znummax = zr
399	ELSE
400	znummax = ( zr + zl ) / 2.
401	zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
402	znum = znummax - 1.
403	zws = xkr_wsbio_min * xkr_zeta / zdiv &
404	& - ( xkr_wsbio_max * xkr_eta * znum * &
405	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
406	& - wmax
407	IF( zws * zwl < 0. ) THEN
408	zr = znummax
409	ELSE
410	zl = znummax
411	ENDIF
412	zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
413	znum = zl - 1.
414	zwl = xkr_wsbio_min * xkr_zeta / zdiv &
415	& - ( xkr_wsbio_max * xkr_eta * znum * &
416	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
417	& - wmax
418
419	zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
420	znum = zr - 1.
421	zwr = xkr_wsbio_min * xkr_zeta / zdiv &
422	& - ( xkr_wsbio_max * xkr_eta * znum * &
423	& xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
424	& - wmax
425
426	IF ( ABS ( zws ) <= xacc ) EXIT iflag
427
428	ENDIF
429
430	END DO iflag
431
432	xnumm(jk) = znummax
433	WRITE(numout,*) ' jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
434
435	END DO
436
437	END SUBROUTINE p4z_sink_init
438
439	#else
440
441	SUBROUTINE p4z_sink ( kt, jnt )
442	!!---------------------------------------------------------------------
443	!! * ROUTINE p4z_sink *
444	!!
445	!! ** Purpose : Compute vertical flux of particulate matter due to
446	!! gravitational sinking
447	!!
448	!! ** Method : - ???
449	!!---------------------------------------------------------------------
450	INTEGER, INTENT(in) :: kt, jnt
451	INTEGER :: ji, jj, jk
452	INTEGER :: iksed
453	REAL(wp) :: zagg1, zagg2, zagg3, zagg4
454	REAL(wp) :: zagg , zaggfe, zaggdoc, zaggdoc2
455	REAL(wp) :: zfact, zstep, zwsmax
456	#if defined key_trc_dia3d
457	REAL(wp) :: zrfact2
458	#endif
459	CHARACTER (len=25) :: charout
460	!!---------------------------------------------------------------------
461
462	zstep = rfact2 / rjjss ! Timestep duration for biology
463
464
465	! Sinking speeds of detritus is increased with depth as shown
466	! by data and from the coagulation theory
467	! -----------------------------------------------------------
468
469	iksed = 10
470
471	DO jk = 1, jpkm1
472	DO jj = 1, jpj
473	DO ji=1,jpi
474	zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
475	wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
476	END DO
477	END DO
478	END DO
479
480	! LIMIT THE VALUES OF THE SINKING SPEEDS
481	! TO AVOID NUMERICAL INSTABILITIES
482
483	wsbio3(:,:,:) = wsbio
484	!
485	! OA Below, this is garbage. the ideal would be to find a time-splitting
486	! OA algorithm that does not increase the computing cost by too much
487	! OA In ROMS, I have included a time-splitting procedure. But it is
488	! OA too expensive as the loop is computed globally. Thus, a small e3t
489	! OA at one place determines the number of subtimesteps globally
490	! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
491
492	DO jk = 1,jpkm1
493	DO jj = 1, jpj
494	DO ji = 1, jpi
495	zwsmax = 0.8 * fse3t(ji,jj,jk) / zstep
496	wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax )
497	wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax )
498	END DO
499	END DO
500	END DO
501
502	wscal(:,:,:) = wsbio4(:,:,:)
503
504	! INITIALIZE TO ZERO ALL THE SINKING ARRAYS
505	! -----------------------------------------
506
507	sinking (:,:,:) = 0.e0
508	sinking2(:,:,:) = 0.e0
509	sinkcal (:,:,:) = 0.e0
510	sinkfer (:,:,:) = 0.e0
511	sinksil (:,:,:) = 0.e0
512	sinkfer2(:,:,:) = 0.e0
513
514	! Compute the sedimentation term using p4zsink2 for all
515	! the sinking particles
516	! -----------------------------------------------------
517
518	CALL p4z_sink2( wsbio3, sinking , jppoc )
519	CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
520	CALL p4z_sink2( wsbio4, sinking2, jpgoc )
521	CALL p4z_sink2( wsbio4, sinkfer2, jpbfe )
522	CALL p4z_sink2( wsbio4, sinksil , jpdsi )
523	CALL p4z_sink2( wscal , sinkcal , jpcal )
524
525	! Exchange between organic matter compartments due to
526	! coagulation/disaggregation
527	! ---------------------------------------------------
528
529	DO jk = 1, jpkm1
530	DO jj = 1, jpj
531	DO ji = 1, jpi
532	zfact = zstep * xdiss(ji,jj,jk)
533
534	! Part I : Coagulation dependent on turbulence
535	! ----------------------------------------------
536
537	# if defined key_off_degrad
538	zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
539	# else
540	zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
541	# endif
542
543	# if defined key_off_degrad
544	zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
545	# else
546	zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
547	# endif
548
549	! Aggregation of small into large particles
550	! Part II : Differential settling
551	! ----------------------------------------------
552
553	# if defined key_off_degrad
554	zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
555	# else
556	zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
557	# endif
558
559	# if defined key_off_degrad
560	zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
561	# else
562	zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
563	# endif
564
565	zagg = zagg1 + zagg2 + zagg3 + zagg4
566	zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
567
568	! Aggregation of DOC to small particles
569	! --------------------------------------
570
571	zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) ) &
572	# if defined key_off_degrad
573	& * facvol(ji,jj,jk) &
574	# endif
575	& * zfact * trn(ji,jj,jk,jpdoc)
576
577	zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc) &
578	# if defined key_off_degrad
579	& * facvol(ji,jj,jk) &
580	# endif
581	& * trn(ji,jj,jk,jpdoc)
582	!
583	! Update the trends
584	! -----------------
585	!
586	tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
587	tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
588	tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
589	tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
590	tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2
591
592	END DO
593	END DO
594	END DO
595
596	# if defined key_trc_dia3d
597	zrfact2 = 1.e3 * rfact2r
598	trc2d(:,:, 5) = sinking (:,:,iksed+1) * zrfact2
599	trc2d(:,:, 6) = sinking2(:,:,iksed+1) * zrfact2
600	trc2d(:,:, 7) = sinkfer (:,:,iksed+1) * zrfact2
601	trc2d(:,:, 8) = sinkfer2(:,:,iksed+1) * zrfact2
602	trc2d(:,:, 9) = sinksil (:,:,iksed+1) * zrfact2
603	trc2d(:,:,10) = sinkcal (:,:,iksed+1) * zrfact2
604	# endif
605	!
606	IF(ln_ctl) THEN ! print mean trends (used for debugging)
607	WRITE(charout, FMT="('sink')")
608	CALL prt_ctl_trc_info(charout)
609	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
610	ENDIF
611
612	END SUBROUTINE p4z_sink
613
614	#endif
615
616	SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra )
617	!!---------------------------------------------------------------------
618	!! * ROUTINE p4z_sink2 *
619	!!
620	!! ** Purpose : Compute the sedimentation terms for the various sinking
621	!! particles. The scheme used to compute the trends is based
622	!! on MUSCL.
623	!!
624	!! ** Method : - this ROUTINE compute not exactly the advection but the
625	!! transport term, i.e. div(u*tra).
626	!!---------------------------------------------------------------------
627	INTEGER , INTENT(in ) :: jp_tra ! tracer index index
628	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed
629	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe
630	!!
631	INTEGER :: ji, jj, jk, jn
632	REAL(wp) :: zigma,zew,zstep,zign, zflx
633	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz
634	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink2
635	!!---------------------------------------------------------------------
636
637	zstep = rfact2 / 2.
638
639	ztraz(:,:,:) = 0.e0
640	zakz (:,:,:) = 0.e0
641
642	DO jk = 1, jpkm1
643	# if defined key_off_degrad
644	zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1) * facvol(:,:,jk)
645	# else
646	zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1)
647	# endif
648	END DO
649	zwsink2(:,:,1) = 0.e0
650
651
652	! Vertical advective flux
653	DO jn = 1, 2
654	! first guess of the slopes interior values
655	DO jk = 2, jpkm1
656	ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
657	END DO
658	ztraz(:,:,1 ) = 0.0
659	ztraz(:,:,jpk) = 0.0
660
661	! slopes
662	DO jk = 2, jpkm1
663	DO jj = 1,jpj
664	DO ji = 1, jpi
665	zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
666	zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
667	END DO
668	END DO
669	END DO
670
671	! Slopes limitation
672	DO jk = 2, jpkm1
673	DO jj = 1, jpj
674	DO ji = 1, jpi
675	zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * &
676	& MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
677	END DO
678	END DO
679	END DO
680
681	! vertical advective flux
682	DO jk = 1, jpkm1
683	DO jj = 1, jpj
684	DO ji = 1, jpi
685	zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
686	zew = zwsink2(ji,jj,jk+1)
687	psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
688	END DO
689	END DO
690	END DO
691	!
692	! Boundary conditions
693	psinkflx(:,:,1 ) = 0.e0
694	psinkflx(:,:,jpk) = 0.e0
695
696	DO jk=1,jpkm1
697	DO jj = 1,jpj
698	DO ji = 1, jpi
699	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
700	trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
701	END DO
702	END DO
703	END DO
704
705	ENDDO
706
707	DO jk=1,jpkm1
708	DO jj = 1,jpj
709	DO ji = 1, jpi
710	zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
711	trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + 2. * zflx
712	END DO
713	END DO
714	END DO
715
716	trn(:,:,:,jp_tra) = trb(:,:,:,jp_tra)
717	psinkflx(:,:,:) = 2. * psinkflx(:,:,:)
718
719	!
720	END SUBROUTINE p4z_sink2
721
722	#else
723	!!======================================================================
724	!! Dummy module : No PISCES bio-model
725	!!======================================================================
726	CONTAINS
727	SUBROUTINE p4z_sink ! Empty routine
728	END SUBROUTINE p4z_sink
729	#endif
730
731	!!======================================================================
732	END MODULE p4zsink

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

Download in other formats: